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PART I 


THE TWO DIMENSIONAL PROBLEM (GRID) USING 
THE FFT-CONJUGATE GRADIENT METHOD 



ABSTRACT 


In some applications, the wires used to construct the 
grids are plated over with highly conducting materials such as 
gold or silver. In those cases, depending on the frequency of 
operation, the coating may not be thick enough to prevent 
currents from flowing in the substrate. The conjugate gradient 
method, in conjunction with the fast Fourier transform is 
employed to solve the problem of scattering from such rectangular 
grids. An internal impedance is utilized to account for the 
effects of the substrate conductivity on the induced current 
densities. Calculated values of the reflection coefficient and 
induced currents for different coating thicknesses, angles of 
incidence and polarizations are presented and discussed. 
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INTRODUCTION 

The subject of reflection and transmission of plane 
electromagnetic waves from grids made of perfectly conducting 
wires, or wires with finite conductivity has been studied by a 
number of investigators [1-5]. In some applications, such as, 
mesh deployable antennas, the mesh wires are made of molybdenum 
substrate coated with a highly conductive material. Depending on 
the frequency of operation, the depth of penetration for the 
incident wave can be larger than the coating thickness. In these 
cases, the electric field will penetrate in the substrate. The 
effects of the resistivity of the substrate on the reflection 
coefficient and induced currents are the aims of this study. 

The FFT-con jugate gradient method [6-12] is used to solve for 
the induced currents on the conducting strips of the grid. An 
internal impedance is used to account for any losses due to the 
finite conductivities of both the substrate material and the 
coating material. This impedance is a function of a) coating 
thickness, b) frequency of operation, c) conductivity of coating 
material, and c) conductivity of substrate material. 

In this work, results for single infinite grids are given as 
a function of the substrate conductivity, coating thickness, and 
polarization of the incident wave. The effects of the substrate 
material on the reflection coefficient and induced currents are 
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also discussed for different angles of incidence. This study can 
easily be extented to the problem of any number of cascaded 
grids. In the case of cascaded grids, there is an interference 
action [13] which is usually extremely sensitive to the losses of 
the wires, which in this particular case it would primarily be 
due to the finite conductivity of the substrate conductivity. 
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2. CURRENT DENSITY FORMULATION 

The magnetic field H due to an electric current density J 
is given by : 


H ( x , y ) = 


V x A ( x , y , z ) 


( 2 . 1 ) 


where A is the associated magnetic vector potential. A and J 
are related by the free space Creen's function 


G(r) = 


A A 

exp(-jk . r) 
4 7T r 


as follows: 


A(?) 


b 



J(r') 


( 2 . 2 ) 


From this the electric field intensity E can be derived from 
Maxwell's equations and expressed as: 

_ VV* T(x,y,z) 

E (x,y,z) = -j w A(x,y,z) + (2.3) 

j co jj, e 

For planar structures we set the z-component of the magnetic 
vector it equal to zero. Now upon expanding equation (2.3) in 
cartesian coordinates we obtain, for z=0: 
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(2.4) 


Considering the periodicity of the two dimensional strucuture 
shown in Figure (2.1) (planar structure), and taking the Fourier 
trasform of equation (2.4) lead to : 

k - a - ot 3 

-t o mn mn H mn ^ 


E ( 


a , 

mn 




jute 


- a 


mn mn 


k 2- 3 
o H mn 


(2.5) 

where the sign (~) denotes the Fourier trasformed quantity. 

a and (3 represent the Floquet coefficients which are 
mn mn 

defined as: 


a: = 2 77 m/a - k sin# sin0 and 

mn o ^ 

= 2 77 n/c -2 m/a cot.Q -k Q sin# sin0 

G( a , (3 ) =- j/2 (k 2 - a 2 -3 2 ) 1/2 7 

v mn H mn o mn P mn 

is the Fourier trasform of Green's function, and J x , are the 

~>s 

unknown current densities .Notice that the spectrum of E . is 
discrete,, that is, it exists for discrete values of a mn and 
<3^. Note , a Iso, that the convolution problem is avoided and instead 
of dealing with an integro differential equation we only 
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have to consider algebraic equations. Taking the inverse Fourier 
trasform of equation (2.5) yields: 


t s (x,y) 


~!_V 


jcoe 


L _ 

mn 


k 2 - 2 

o a mn 


- a. B 

mn p mn 


a mn ^mn 


. 2 2 
k o " 0 


mn 


G. J 


■exp [ j ( a x+ By) ] 
mn mn 


( 2 . 6 ) 


To enforce the boundary condition over the surface of all 
metallic regions we require that the total tangential electric 
field should satisfy the condition : 


E 1 (x,y ) +'E s (x,y) =0 


(2.7) 


where E is the incident electric field and E is the scattered 

• • — . 
electric field. Substituting for the value of E from equation 

(2.7) into equation (2.6) yields: 



fl mn ^mn 



G(a__ / B mr , ) J (a mn /|3mn ) 
mn ^mn mn f-'mn 

.exp [ j ( a x+ B mn y ) 

mn ^mn (2>8) 


Equation (2.8) can be recognized as the inverse discrete 
Fourier transform which can be performed via the fast Fourier 



transform (FFT). Equation (2.8) could now be written in an 
operator form as: 
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where Z is the product of G, the Floquet modes and the inverse 
Fourier transform. 

A solution of the above equation will yield the unknown 
current densities J x and from which the reflected and 
transmitted fields can be obtained and hence the reflection and 
transmission coefficients could be calculated. 

Now, one way to solve for and is to use the conjugate 
gradient method [12]. To guarantee a convergent scheme, equation 

(2.9) has to be properly modified. To do that, multiply both 

★ 

sides of equation (2.9) by Z (i.e. the conjugate transpose of 

Z ) to obtain: 
mn 


★ — * — y 

-Z E = Z Z J 
mn mn mn 


( 2 . 10 ) 


★ 

where the product Z Z is a Hermitian matrix and therefore 
^ mn mn 

positive definite. Now the conjugate gradient method can be 


applied directly to equation (2.10). 
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3. FORMULATION OF THE FFT-CON JUGATE GRADIENT METHOD FOR 

COATED WIRES 


3.1 EQUIVALENT RADIUS PRINCIPLE AND INTERNAL IMPEDANCE 


The strip analysis can be used to determine the scattering 
characteristics from a mesh of cylindrical wires by employing the 
"equivalent radius principle" concept. This is accomplished by 
replacing the non-circular cross section of a metallic strip with 
a circular wire whose radius is the "equivalent radius" of the 
non-circular cross section as shown in Figure (3.1). Butler [14] 
has shown that the equivalent radius of a narrow conducting strip 
is one fourth of its width i.e 

a^„ =a/4 
eq 

where a^ is the equivalent radius of a cylindrical wire, and a 
is the width of a thin metallic strip. 



Fig. 3.1 Equivalent radius of a strip 



- 8 - 

For the case where the wires are coated (See Figure 3.2) 
the necessary boundary condition that must be satisfied is : 


E 


i 



(3.1) 


instead of E +E =0 , where I is the current in the wires and 

Zint is the internal impedance of the wire. For coated wires 

Z. . is given by [15] : 
int 


(1+j) (sinh(t 1 d) + ( r s 2 / r s ].) cosh(tjd)] 


b [cosMt^d) + ( R s 2/ R sl^ sinh(t-^d)] 


where t n = ( 1+j ) / nf JU ^ 


and 


(3.3) 


j ) 


t,“(l+j ) TTf A 2 °2 


(3.4) 


d=coating thickness 


R sl = V 7Tf CT 1 and R s2 = V ?Tf W 


(3.5) 


/i^= permeability of coating 
cr^= conductivity of. coating 
pi 2 = permeability of substrate 


0 2 = conductivity of substrate 




e 
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This expression for z ^ nt can now be used in equation (2.9)i.e, 


*— i — > — 

E + ET = Z . . I = (Z. A ) J = Z. ^ J 

1 f- ' 1 4- ' i m 4- 


int 


mt 


int 


(3.6) 


since J=I/A where A is the surface area of the wire. This leads 
to : 



(3.7) 


Replacing this expression for E in equation (2.6) yields : 

s 


1 , ^ 


-E + Z 


int 


— > 

J = 


mn 


or 



(Z 

mn 



) J 


B J 


(3.8) 


Now equation (3.8) can be solved for J. using the conjugate 

gradient algorithm. Rather than form the matrix (Z -Z . . ) 

mn int 

explicitly, one can apply the conjugate gradient method using the 
algorithm given on the next page. 

In Figures 3.3 and 3.4, and internal impedance with a ratio 
of r s2 / r s i = ®* and 1*6 are shown, respectively [15]. These 

figures are important, because they show that z ^ nt will not be 
that of the coating material alone, despite the fact that the 
ratio of the coating thickness to the skin depth is greater than 


one . 



CONJUGATE GRADIENT ALGORITHM 



The equations for the n iteration are 
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( 3 . 9 ) 



-Hn+l) . z * f(n+l) _ T * — ( n+1 ) + ,j ■*(»> 

H mn int n 


END OF DO LOOP 




Fig. 3.3 Normalized Impedance. Solder coating on a copper substrate. 

Rs2/Rsl =0.34 
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In the current algorithm, if this ratio is less than the 
number four, equation (3.2) is calculated. On the other hand, 
if the ratio is greater than the number four, ^ is due to 

the coating material only, and it is given by [16] : 


m 


int 


2 n a 


eq 


I (7a ) 

o ' eq 




3. 10) 


where 


m 


j \ 

' O + j 

m m 


1/2 


is the intrinsic impedance of the metal. 7 is equal to 

1/2 . . 

(j h co(g + j to e )) and I and I n are the modified 

“m m m o 1 

Bessel functions. At high frequencies the ratio I /Ij. is 
approximately equal to one. Thus, any uneccessary computations 
are avoided by evaluating equation (3.10) instead of equations 
(3.2), (3.3), (3.4), and (3.5). 


3.2 SOLUTION OF APERTURE FIELDS 


To solve for the aperture fields (See Figure 3.5), the 
incident H field is expressed in terms of the electric field E in 
the aperture region by : 
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mn 




'mn ^mn 


G( ft. a __ ) E i 

mn h mn 


• exp [ j ( a m x+ flv) ] 
mn r mn 


( 3 . 11 ) 


h* — BB — H 



Fig. 3.5 Sampling for the Aperture fields. 


For a complete derivation of equation (3.11) see [17]. 
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4. REFLECTION COEFFICIENTS 


The transmission and reflection coefficients are the 
quantities of most importance in characterizing the properties of 
a mesh. In order to define those coefficients for both 
polarizations, transverse electric (TE) and transverse magnetic 
(TM), it is necessary to first define the incident and scattered 
fields. For TE polarization, the incident fields are : 

E = E sin (-0 ) ; E = E cos 0 
x o ' r y o v 



E cos 0 
o 

V 


COS 0 



E sin0 
o _ 

r? 


cos 0 


where E is the amplitude of the incident electric field and 
o 

1/2 

77 -{floj € o ) ' is the free space wave impedance. For TM 
polarization, the incident fields are given by : 


E x = E q COS0 COS0 

E sin ( 0 - 7T /2 ) 
H - — 2 

x 1 


E 

y 


=E q cos 9 sin0 
E q cos ( 0 -71/2 ) 
7? 


According to Wait and Hill [ 4 ] , when the spacing between 

adjacent wires of the mesh is less than A/2, only the J and 

J components contibute to the scattered field. and 

ooy c oox 

J ooy are the zero-mo< 3e current density components. The 

rectangular components of the scattered field, can be obtained 
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from equation (2.6) as follows : Solve for J and substitute the 


solution in equation (2.6) to obtain the E and E components 

X y * 

of the scattered field. Once these components are found the 
reflection (amplitude) coefficient becomes : 

(4.1) 

(4.2) 

If the total power reflection coefficient R is desired then the 
following expression can be used : 

Real 


R = E S /(E 2 + E 2 ) 1//2 

x x ' x y 

R = E S /(E 2 + E 2 ) 1/2 

y y ' x y 


R 


“e S x H S 


A 

z ds 


'-unit cell 


REAL 


E 1 x H 1 


unit cell 


(-z) ds 


(4.3) 


H s is the scattered magnetic field derived from E by making use 
of Maxwell's equations. 

Moreover, if the total power transmission coefficient, T, is to 
be computed, one can employ the formula below : 

Real , 

, aperture 


E a x H a . (-z) dA 


(4.4) 


Real J JE 1 x H 1 . (-z) dA 
l aperture 

where E a is the aperture electric field and H a is the aperture 
magnetic field. 
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5. RESULTS 

A number of cases were examined to check the reflection 
coefficients and induced current densities as functions of wave 
polarization, angle of incidenece, wire spacing, coating 
thickness, and substrate conductivity. Since no other 
theoretical or experimental data currently exist, only results 
obtained via this new algorithm are presented herein. 

Figure (5.1) depicts the change in the reflection 

coefficient (amplitude) as the thickness of the coating material 

changes. In this case, the substrate conductivity is 50 S/m and 

0 

the coating conductivity is 10 S/m. As expected, the thinner the 
coating material is the deeper the fields will penetrate and 
hence more losses should be expected in the strips. That means 
that the amplitude of the reflection coefficient will decrease. 
This figure also shows that when the thickness of the coating 
material is large the losses are mainly due to the finite 
conductivity of this coating material. On the other hand, at 
small coating thicknesses the losses will be primarily due to the 
substrate conductivity. Figure (5.2) depicts the behavior of the 
reflection coefficient for two different angles of incidence as 
the substrate conductivity is varied. The thickness and 
conductivity of the coating material are kept constant. From 
this figure one can see that as the conductivity of the substrate 
gets larger the reflection coefficient increases. This is an 
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anticipated result since the losses in the strips decrease by 
incresing the conductivity of the substrate. 

Figure (5.3) exhibits a similar behavior in the case of 
thansverse magnetic polarization for angles 0 = 60° and 0=0°. By 
comparing Figures (5.2) and (5.3), one can make the additional 
observation that for TM polarization the reflection coefficient 
does not change as drastically as in the TE polarization case. 

5 

Table (5.1) shows that if the substrate conductivity is 10 S/m 

0 

and the coating conductivity 10 S/m , the change in the 
reflection coefficient is not very much. The fact that the 
reflection coefficient for 0=0° is larger than the reflection 

O 

coefficient for 0=70 is due to the general behavior of 
conducting grids ( a = , or finite) and not due to the coated 
material. In Table (5.2) the reflection coefficient is calculated 

_ O O 

for TM polarization at 0=0 and 0=70. Note that the substrate 
conductivity is much smaller than the conductivity of the the 
coating material .This difference is basically the reason for 
having a noticable change in the reflection coefficient at small 
coating thicknesses. 

Next, the current densities were calculated to study their 
behavior as we vary the substrate conductivity and the angle of 
incidence . Figures (5.4) and (5.5) depict the behavior of 

the copolar component of the current density J^, for a grid with 
wide strips. The current densities do not really change very 
drastically when the substrate conductivity is changed. The 
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e 





FIG. 5.3 REFLECTION COEFFICIENT VERSUS SUBSTRATE CONDUCTIVITY 
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TABLE 5.1. Reflection coefficients for a cell of a=0.25\ 


and a strip 

width 

of 0.005^ , with 9 ■ 

=0 and 

QS 

II 

-J 

o 

o 

g 

respectively. a =10 

* 1 coat 

S/m and 

o K =io 5 

sub 

S/m . 

TE polarization case 

with a 

sampling rate of 

32x32 samples 


Thickness 


Reflection 

Reflection 

of coating 


coefficient 

o 

coefficient 
. o 

in X 


0 =0 

9 =70 

10- 3 


0.6073501 

0. 9077239 

io' 4 


0. 6073501 

0. 9077239 

io“ 5 


0.6086803 

0.9076653 

10- 6 


0. 6066462 

0. 9077298 

io" 7 


0.6063545 

0. 9064933 

00 

l 

o 
« — 1 


0. 6060203 

0. 9045778 

io' 9 


0. 6028293 

0.9045168 

H)' 10 


0.6022630 

0.9045149 
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Table 5.2. Reflection coefficients for TM incidence and a 

square cell of a=0.25^.The strip width is 0.005^ 

O =10® S/m and a , =10 3 S/m. The 
coat sub / 

sampling rate is 16x16 samples. 


Thickness 

Reflection 

Reflection 

of coating 

coefficient 

o 

coefficient 
„ o 

in \ 

9 =0 

0=70 

10 -3 

0. 7405149 

0. 4434715 

i 

o 

i — t 

0.7405149 

0.4434715 

10 -5 

0.7405282 

0. 4434533 

10- 6 

0.7405004 

0.4434442 

10 -7 

0. 7385695 

0. 4423161 

H-I 

o 

i 

00 

0.7225416 

0.4307751 

10~ 9 

0.6665412 

0. 3766584 

o 

1 — t 

1 

o 
1 — 1 

0.6449590 

0.3511487 




FIG. 5.3 CURRENT DENSITIES 
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general behavior, though, of the current density on the 
strip, i.e. J is large at the edges and it dips at the center, 
is still maintained. 

Finally, the reflection coefficients of the structure shown in 
Figure (5.6) were studied , since this structure generates a 
cross polarized component. Tables (5.3) and (5.4) show that for 
both polarizations, TM and TE, the changes in the thickness of 
the coating material do not correspond to any drastic changes in 
the amplitude of the reflection coefficient. This is an 
interesting observation since the actual weaved mesh used for 
deployable antennas resembles this last structure more than the 
rectangular periodic structure. In both of these tables the wave 
was normally incident on the periodic structure . The substrate 

3 

conductivity was a =10 S/m and the conductivity of the coating 

9 

material was o =10 S/m. 




- 29 - 


TABLE 5.3 Reflection coefficients for a cell of a=b=0.55 


and a strip width of 0.005\ with 0=70 and 

o 9 

0=80 , respectively. a CO at = ^ s / m ' and 

3 

<7 s ub = 10 S/m. A TM polarization case with 
a sampling rate of 32x32 samples is used. 


Thickness 

Reflection 

Reflection 

of coating 

coefficient 

coefficient 

in > 

0 

o 

o 

r-' 

n 

0 

o 

= 80 


copolar 

crosspolar 

copolar 

crosspo 

10~ 5 

0.3372 

0. 1074 

0. 3169 

0.0517 

10~ 7 

0.3370 

0.1074 

0.3169 

0.0517 

CO 

1 

o 

f— I 

0.3368 

0. 1072 

0. 3166 

0.0516 

10~ 9 

0.3338 

0.1051 

0.3136 

0.0508 

10' 10 

0.3273 

0. 1023 

0.3073 

0.0498 
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TABLE 5.4 Reflection coefficients for a cell of a=b=0.55 

and a strip width of 0.005)^ , with 0=70° and 

_ ' o q 

Q =80 , respectively. O coat = 10 s / m and 

3 

a sub =1 ° s / m * A TE polarization and a sampling 
rate of 32x32 samples is used. 

Thickness Reflection Reflection 


of 

coating 

coefficient 

coefficient 

in 

X 

e 

= 70° 

9 =80° 




copolar 

crosspolar 

copolar 

crosspolar 

10 -5 


0.28756 

0. 10101 

0.30488 

0.05212 

10" 7 


0. 28755 

0. 10103 

0.30493 

0.05235 

00 

1 

o 

r~H 


0. 28726 

0. 10083 

0.30462 

0.05210 

10" 9 


0. 28441 

0.09671 

0.30171 

0.05151 

10- 10 


0.27840 

0.09805 

0.29547 

0.05064 
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6. CONCLUSIONS AND RECOMMENDATIONS 


An algorithm was developed which can solve the problem of 
electromagnetic scattering from meshes made of coated wires. The 
algorithm solves for the reflection coefficients and the induced 
current densities. A variety of grids were checked as functions 
of polarization, angle of incidence , coating thickness coating 
conductivity, and substrate conductivity. Generally, it was found 
that when the difference in conductivity between the coating 
material and the substrate material is large the algorithm is 
very sensitive to the losses due to the finite substrate 
material. Also, when the thickness of the coating material 
becomes small compared to the depth of penetration the losses are 
primarily due to the substrate material. On the other hand, if 
this thickness is comparable or larger than the depth of 
penetration the losses are determined by the finite conductivity 
of the coating material. Moreover, it was observed that for non 
rectangular structures the amplitude of the reflection material 
does not change very much even at small coating thicknesses. 

This algorithm could be extended to more complicated 
structures such as the skew-symmetric grids and cascaded grids. 
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8. LISTING OF THE FFT-CG METHOD FOR THIN STRIPS 


C ****C ON JUGATE GRADIENT METHOD ***** 

C **** SOLVES FOR THE CURRENT DENSITIES *** 

C **** MINIMIZATION IN THE RANGE **** 

COMPLEX CONE. CZERO. CXMN, CREFX. CREFY. CRET. CREF. ZINT. GUESS. Tl. 
. SINHH, COSHH. CTT , CTF 
COMPLEX YC32, 32>/1024*(0. 0, 0. 0)/ 

COMPLEX X < 32, 32 ) /1024* (0. 0.0. 0 ) / 

COMPLEX G < 32, 32 ) /1024* (0. 0,0. 0 ) / 

COMPLEX XU<32, 32)/1024*<0. 0, 0. 0)/ 

COMPLEX YU<32, 32)/I024*(0. 0, 0. 0)/ 

COMPLEX RX(32, 32)/1024*<0. 0, 0. O)/ 

COMPLEX RY ( 32, 32 ) /1024* (0. 0. 0. 0 ) / 

COMPLEX J, HXI, HYI,CWK<32>, F10 
COMPLEX DY<32, 32)/1024*<0. 0, 0. 0)/ 

COMPLEX DX ( 32, 32 > /1 024* ( O. 0, 0. 0 ) / 

COMPLEX T X ( 32 , 32)/1024*<0. 0, 0. 0)/ 

COMPLEX TY<32, 32)/1024*(0. 0. 0. O)/ 

REAL K, K2, RWK(342> 

DIMENSION IWK< 342), RR<250>, CH<250 > , AMP < 32 > . RINDEX < 32) , 

. CRDS < 32 ) 

REAL U(32)/32*0. 0/ 

REAL V<32, 32) /I 024*0. 0/ 

OPEN (10. FILE= ' FASTDATA ', STATUS= 'OLD ' > 
open (11, f i le= ' fast out ', s tat us= 'new ' ) 

WR ITE ( *, * ) 'INPUT AA, BB, CC, DD, F, ERR ' 

READOO, *> AA, BB.CC.DD 
C22 FORMAT (8E10. 4) 

F=2. 998E+8 

C **** I0PT=0 FOR A RECTANGULAR MESH ******* 

C **** I0PT=1 FOR A PARALLEL GRID ********** 

C 

IX =32 
I0PT=0 

IF ( IOPT. GT. O) CC=1. 500E+15 
IF ( IOPT. GT. O) DD=1. 500E+15 
WRITE <11, *) AA, BB,CC,DD 
C . 'DD= ', FI 5. 8, 'ERR= '.F15.8) 

WR ITE ( *, * ) F 

C44 FORMAT ('O'. * FREQ = '.E10, 4) 

WR ITE (11. *) 'INPUT PHI. THI, PSI ' 

READ(10, *) PHI, THI, PSI 
WR ITE <11, * ) PHI, THI, PSI 

C55 FORMAT ( 'O'. ' PHI= '.F10. 1. ' THETA= ',F10. 1,' PSI= ',F10. 1) 

C *** ITM=1 FOR TM POLARIZATION ****** 

READ <10, *) ITM 

WR ITE <11, * ) 'INPUT ITM=0 FOR TE POLAR. OR ITM=1 FOR TM POL' 
WR ITE <11, * ) ITM 

WR ITE <11. * ) 'INPUT CONDUCTIVITY OF COATING' 

READUO,*) SIGMA1 
WRITE! 11, *) SIGMA1 

WR ITE <11, * ) 'INPUT CONDUCTIVITY OF SUBSTRATE' 

READUO, *) SIGMA2 
WRITE (11, * ) SIGMA2 

WR ITE (11, * ) 'THICKNESS OF COATING ?' 

READUO,*) THICKNESS 
WRITE (11,*) THICKNESS 
C *** READ THE NUMBER OF ITERATIONS **** 

READUO,*) NO I 

WRITE (11,*) 'NUMBER OF ITERATIONS' 

WR ITE <11, * > NOI 



Cb6 FORMAT! 13) 

PI =3. 141593 
PI2=PI/2. 

TP 1=6. 283185 
CV=2. 997956E+8 
UU=4. E-7*PI 
RTD=57. 29578 
EP=8. 854E-12 
ETA=SQRT ( UU/EP > 

J=CMPLXCO. O. 1. 0 ) 

ITER = 0 

CONE=CMPLX< 1. 0, O. 0) 

CZERO=CMPLX (0. O, 0. 0) 

W=TPI*F 

C *** COMPUTE INTERNAL IMPEDANCE OF STRIP **** 

C 

Tl = < 1. O, 1. 0>*SQRTCPI*F*UU*SIGMA1 ) 

RS1=SQRTCPI*F*UU/SIGMA1 ) 

RS2=SGR T ( P I *F*U U/SI GMA2 ) 

SKIND1=1. /SORT! PI*F*UU*SIGMA1 ) 

RATI 0=THICKNESS /SKI ND1 
IF (RATIO. GE. 4. O ) THEN 

ZINT=( 1. O. 1. O )*SGRTCPI*F+UU/SIGMA1 ) 

ELSE 

Z I NT = < 1 . O, 1. O)* <SINHHCT1*THICKNESS)+CRS2/RS1 )*C0SHHCT1*THICKNESS> ) 
Z INT=Z INT *RS1 / < C0SHHCT1*THICKNESS ) + < RS2/RS1 )*SINHHCT1*THICKNESS) ) 
END IF 

WRITE!*, *) RATIO, ZINT, RSI, RS2, Tl, THICKNESS, SKIND1 

ALAMB=CV/F 

AA=AA/ ALAMB 

BB=BB/ALAMB 

CC=CC/ ALAMB 

DD=DD/ ALAMB 

C *** DETERMINE SAMPLING POINTS THAT CORRESPOND TO THE 
C CONDUCTING REGIONS AND THE APERATURE ********* 

NX=IFIXCBB/AA*FL0ATCIX)*2. ) /4*2 
NY=IFI X C DD/CC*FLOAT ( I X ) *2. )/4*2 
NX1--(IX-NX)/P+1 
NX2=NX1+NX-1 
NYl=CIX-NY)/2 +1 
NY2=NY1+NY-1 

WRITE <11, *) NX, NX 1, NX2, NY, NY1 , NY2 

K=TP I /ALAMB 

K2=K**2 

STSPK=SIN <THI/RTD)*SINCPHI/RTD)*K 
STCPK=S I N ( THI /R TD ) *COS ( PHI /RTD ) *K 
CPS=COS(PSI/RTD ) /SIN ( PSI /RTD ) 

70 CONTINUE 

C *** CALCULATE FLOQUET MODES ****** 

DO 100 M=l, IX 
IF <M. GT. IX/2+1 ) GOTO 75 
U < M ) =TP I * < M-l ) / AA-STCPK 
GOTO 80 

75 U CM) = TP I* (M— I X-l ) /AA-STCPK 

80 CONTINUE 

DO 90 N=l. IX 

IF (M. GT. IX/2+1. AND. N. GT. IX/2+1 ) GOTO 84 
IF CM. GT. IX/2+1) GO TO 83 
IFCN. GT. IX/2+1 > GOTO SI 

VC M, N )=TP I* CN-l ) /CC—TP I * < M— 1 ) /AA*CPS-STSPK 
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GOTO 85 

81 V(M, N)=TPI*(N-I X-l )/CC-TPI*(M-l )/AA^CPS-STSPK 
GOTO 85 

83 V ( M> N ) =TP I * CN— 1 ) /CC-TP I# (M- 1 X-l > /AA*CPS-STSPK 
GOTO 85 

84 V(M, N)=TP I*(N-I X-l ) /CC-TP I ♦ ( M-I X- 1 ) /AA*CPS— STSPK 

85 IF(K2. GE. U(M)**2+V(M, N)**2) G (M. N ) =-J*SQRT ( K2— ( U(M) -»*2 

♦ V(M, N>*^2 ) ) 

IF (K2. LT. U(M)**2+V<M. N)*#2) G(M, N ) =-SQRT ( U ( M > ♦♦2+V ( M, N)**2-K2> 

♦ CONE 

90 CONTINUE 
100 CONTINUE 

IF ( ITM. GT. 0) GOTO 110 

C ♦ ♦♦ INCIDENT FIELDS FOR TE POLARIZATION ♦♦♦♦ 

EXI = S I N ( — PHI / RTD) 

EYI = COS ( PHI /R TD ) 

HX I=COS (PHI /RTD ) ♦COS (THI /RTD) /ETA 
HYI=SIN ( PHI /RTD ) ♦COS ( THI /RTD ) /ETA 
EF=1 . O 
GOTO 120 

C ♦♦♦ INCIDENT FIELDS FOR Tfl POLORIZATION ♦♦♦♦ 

110 EX I=COS (PHI /RTD MCOS(THI/RTD) 

EYI=SIN(-PHI/RTD)^COS(THI/RTD) 

HYI=SIN (PHI /RTD— P 12 ) /ETA 
HXI=COS (PHI /RTD— P 12 ) /ETA 
ET=1. 0*C0S < THI / RTD) 

120 CONTINUE 

C**^ CALCULATE THE RESIDUAL VECTORS RX AND RY ♦♦♦♦ 

C 

333 DO 200 M=l, IX 
DO 190 N=l. IX 
RX (Mi N)=EXI+X(M, N) 

RY (Mi N)=EYI+Y(M. N) 

IF (M. GE. NX 1 . AND. M. LE. NX2. AND. N. GE. NY1 . AND. N. LE. NY2 ) 

. RX(M. N)=CZERO 

IF (M. GE. NX 1 . AND. M. LE. NX2. AND. N. GE. NY1. AND. N. LE. NY2 ) 

. RY(M. N)=CZERO 

ERROR=ERROR+RX ( M. N)-*CONJG(RX (M, N) >+RY(M, N> ♦CONJG (RY<M, N) ) 

IF (ITER. EQ. 0) F5=F5+RX(M. N) ♦CONUG (RX ( M, N> )+RY(M, N ) ♦CONJG (RY(M, N) ) 
DX (Mi N ) =R X ( M, N ) 

DY (Mi N ) =R Y ( Mi N ) 

190 CONTINUE 
200 CONTINUE 
C 

MULTIPLY THE RESIDUALS BY THE CONJG. TRANS. OF Z TO 
C FIND THE DIRECTION VECTORS DX AND DY *♦♦♦♦♦* 

C ♦♦♦♦ FIND THE FOURIER TRANSFORM OF THE RESIDUALS ♦♦*♦♦ 

C 

CALL FFT3D(DX, I X, IX. IX, IX, 1, 69, IWK, RWK, CWK) 

CALL FFT3D(DY, I X, IX, IX. IX, 1, 69, IWK, RWK, CWK) 

DO 220 M=l. IX 
DO 210 N=l. IX 
CXMN=DX(M,N) 

DX (M. N ) = ( CONJG ( G(M. N)-V(M. N)^2/G(M, N) ) ♦DX ( M, N>- 
. CONJG ( V( M, N ) *U ( M )/G (M, N) >^DY(M, N) > /CONJG (J*W*EP ) /2. 

DY(M, N ) <= ( CONJG ( -V(M, N)^U(M)/G(M, N) )^CXMN+ 

. CONJG ( G ( M, N)— U (M)**2/G(M> N ) ) ♦DY ( M, N) ) /CONJG ( J^W*EP > /2. 

210 CONTINUE 
220 CONTINUE 
C 
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CALL FFT3D(DX, I X. IX. IX. IX. 1. -69, IWK, RWK, CWK) 

CALL FFT3D ( DY, I X, IX. IX, IX. 1, -69. IWK, RWK, CWK) 

C*** STORE DX AND DY IN TX AND TY ***«• 

DO 360 (1=1. IX 
DO 350 N=l. IX 

DX (M, N)=DX (M, N) -RXCM, N)*CONJG < ZINT) 

DY ((1, N ) =DY ( (1, N) — RY<M, N > *CON JG ( Z I NT ) 

IF CM. GE. NX 1 . AND. M. LE. NX2. AND. N. GE. NY 1 . AND. N. LE. NY2) 

. DX(M, N)=CZERO 

IF <M. GE. NX 1 . AND . M. LE. NX2. AND. N. GE. NY1 . AND. N. LE. NY2) 

. DY(M, N)=CZERO 
TY (M, N)=DY(M, N) 

TX <M, N)=DX (M, N) 

F3=F3+C0NJG(DX( (1, N) )*DX(M, N) +CONJG ( DY ( (1, N) )*DY(M, N) 

350 CONTINUE 
360 CONTINUE 

C***« THE ITERATIVE PROCESS STARTS NOW ! ! ! ! ! *** 

365 CALL FFT3D(TX, I X, IX, IX, IX, 1, 69, IWK, RWK, CWK) 

CALL FFT3DCTY, I X, IX, IX, IX, 1, 69, IWK, RWK, CWK) 

DO 400 M=l, IX 
DO 370 N=l. IX 
CXMN=TX (M, N) 

TX(M, N) = ( (G(M, N ) — V( M, N)**2/G(M, N) >*TX(M, N)-(U<M)*V(M, N)/G(M, N) ) 
. *TY(M, N) )/(J*W*EP)/2. 

TY(M, N) = (-U(M)*VCM, N)/G<M, N ) *CXMN+ < G ( M, N)-UCM)**2/G(M, N) ) 

. *TY(M, N) )/<JwW*EP)/2. 

370 CONTINUE 
400 CONTINUE 

CALL FFT3D ( TX, I X, IX, IX. IX, 1,-69, IWK, RWK, CWK) 

CALL FFT3DCTY, I X, IX, IX, IX, 1, -69, IWK, RWK, CWK) 

F1=0. O 

DO 410 M=l, IX 
DO 410 N=l. IX 

TX(M. N)=TX(M, N ) — DX< M, N)*ZINT 
TY (M. N)=TY(M, N) -DY(M, N)*ZINT 

IF (M. GE. NX1 . AND. M. LE. NX2. AND. N. GE. NY1 . AND. N. LE. NY2) 

. TX(M, N) = CZERO 

IF(M. GE NX 1 . AND M. LE. NX2. AND. N. GE. NY1 . AND. N. LE. NY2) 

. TY(M. N)=CZERO 

410 F1=F1+C0NUGCTX ( M, N) )*TX <M, N ) +CONJG ( TY ( (1, N) )*TY<M, N) 

C *** CALCULATE THE FACTOR AN *** 

ITER=ITER+1 

AN=F3/F1 

CH< ITER )=SQRT(abs (error) >/SQRT(F5) 

C *** CALCULATE THE ERROR *** 

ERR0R=ERR0R-(F3**2/F1 ) 

C *** UPDATE THE VALUES FOR X AND Y *** 

DO 560 (1=1- IX 
DO 550 N= 1 . IX 
X<M> N)=X(I1. N)+AN*DX(M, N) 

Y < M, N)=Y(f1. N)t-AN*DY(M, N) 

550 CONTINUE 
560 CONTINUE 

C **** FIND A NEW ESTIMATE FOR THE RESIDUAL VECTORS RX AND RY *** 

DO 5B0 M=l, IX 
DO 570 N=l, IX 

RX(M, N)=RX(M, N) -AN*TX(M, N) 

RY(M. N)=RY(M, N) -AN*TY(M. N) 

TX (M, N)=RX(M, N) 

TY<M, N)=RY(M. N) 



570 CONTINUE 
580 CONTINUE 

RR < I TER ) =FLOAT < I TER > 

WR ITE < *, * ) CH< I TER) , RR < ITER > 

C****#MULT I PLY TX AND TY BY THE CONJG. TRANS. OF THE MATRIX Z ** 
CALL FFT3D ( TX. I X, IX, IX, IX, 1, 69, IWK, RWK, CWK) 

CALL FFT3D < TY, I X, IX, IX, IX, 1, 69, IWK, RWK, CWK) 

DO 600 M=l. IX 
DO 590 N=l, IX 
CXMN=TX<M, N> 

TX (M, N ) = ( CONJG ( G ( M, N)-V(M, N)**2/G<M, N) )*TX(M, N>- 
. CONJG ( V ( M, N)*U <M)/G<M, N) >*TY(M, N) ) /CONJG ( J*W*EP ) /2. 

TY (M, N) = ( CONJG ( -V(M, N ) *U < M ) /G C M, N ) >*CXMN+ 

. CONJG ( G ( M, N ) — U <M)**2/G(M, N) >*TY(M, N) )/CONJG< J*W*EP >/2. 

590 CONTINUE 
600 CONTINUE 

CALL FFT3D ( TX, I X, IX, IX, IX. 1,-69, IWK, RWK, CWK) 

CALL FFT3D ( TY, I X, IX, IX, IX. 1, -69, IWK, RWK, CWK) 

F2=F3 

C F3=F3/400. 

c IF(CH( ITER). LT. 0. 30) F3=0. O 

F3=0. O 

DO 644 M=1 > IX 
DO 644 N-X. IX 

TX (M, N)=TX (M, N ) — RX ( M, N ) *C ON JG ( Z I NT ) 

TY(M, N ) =TY < M, N) — RY(M, N ) *CONJG ( Z INT ) 

IF <M. GE. NX1. AND. M. LE. NX2. AND. N. GE. NY1. AND. N. LE. NY2) 

. TX<M, N)=CZERO 

IF <M. GE. NX 1 . AND. M. LE. NX2. AND. N. GE. NY1. AND. N. LE. NY2 ) 

. TY(M. N)=CZERO 

F3=F3+C0NJG<TX< M, N) )*TX(M, N)+CONJG(TY(M, N) )*TY(M, N) 

644 CONTINUE 

C **** CALCULATE THE FACTOR BN **#* 

BN= < F3/F2 ) 

C ***** UPDATE THE DIRECTION VECTORS DX AND DY ***** 

C 

DO 664 M=l. IX 
DO 654 N=l, IX 

DX(M, N)«TX(M, N) +BN*DX(M, N) 

DY (M, N)«TY(M, N) +BN*DY(M, N) 

TX(M. N)=DX(M, N) 

TY (M. N)=DY(M, N) 

654 CONTINUE 
664 CONTINUE 

C **** GO FOR ANOTHER ITERATION IF YOU WANT **** 

IF ( ITER. GT. NDI ) GO TO BOO 
GO TO 365 

800 DO 820 1=1. IX 

AMP ( I )=CA3S<Y < I , 16) ) 

R INDEX ( I ) = (FLOAT ( I-I X/2 ) — . 5)/IX*AA*l. 045 
WR ITE ( *, * ) AMP ( I ) , R 1 NDE X ( I ) 

820 CONTINUE 

WRITE (11, * ) ITER 

CALL FFT3D ( X, IX, IX- IX, IX, 1, 69, IWK, RWK, CWK) 

CALL FFT3DCY. I X , I X, I X, I X, 1 , 69, IWK, RWK, CWK) 

DO 840 M=l. IX 
DO 830 N= 1 , IX 
CXMN=X(M, N) 

X(M, N) = ( (G(M, N) -V(M. N)**2/G(M, N) )*X (M, N)-(U<M)*V(M. N)/G(M, N) ) 
. *Y <M, N) )/ ( J*W*EP)/2. 



YCM, N)=(-U<M)*V in, N)/G(M, N } « C XMN+ ( G < M, N > -U ( M ) **2/G < M. N> ) 
. *Y(M, N) )/<J*W*EP)/2. 

030 CONTINUE 
040 CONTINUE 

C **** CALCULATE REFLECTION COEFFICIENTS 
CREFX=X ( 1 . 1 ) /FLOAT ( I X ) *#2 
CREFY=Y< 1 , 1 ) / FLOAT ( I X ) **«2 

CREF=CREFX*SIN< -PHI /RTD > +CREFY*C0S < PHI /RTD ) 
CRET=CREFX*COS( PHI/RTD)+CREFY*SIN<PHI/RTD> 

IF( ITM. GT. 0> GO TO 850 
C 

C **** TE POLARIZATION ***• 

C** THIS IS THE CO-POLARIZED COMPONENT **** 

REFF=CABS (CREF/EF) 

CTF=1— CREF 
TFF=CABS( CTF) 

C *** THIS IS THE CROSS-POLARIZED COMPONENT *** 

REFT=CABS ( CRET/ EF ) 

WRITE (II# *) 'CREF. CRET, REFF, REFT » TFF ' 

WRITE < 11.*) CREF, CRET. REFF, REFT. TFF 
050 CONTINUE 
C *** TM POLARIZATION *** 

C *** THE CO-POLARIZED COMPONENT ** 

RETT=CABS ( CRET/ ET ) 

CTT=1— CRET 
TTT=CABS ( CTT) 

C **** THE CROSS-POLARIZED COMPONENT*** 

RETF=CABS < CREF/ ET ) 

WRITEU1, *) 'CRET, CREF. RETT, RETF, TTT ' 

WRITE( 11, *) CRET, CREF, RETT, RETF, TTT 
900 STOP 
END 
C 
C 

COMPLEX FUNCTION SINHH(X) 

COMPLEX X 

SINHH=0. 5* ( CEXP (X)-CEXP(-X ) ) 

END 

C 

c 

COMPLEX FUNCTION COSHH(X) 

COMPLEX X 

C0SHH=0. 5* ( CEXP < X )+CEXP (-X ) ) 

END 



PART II 


THE ONE DIMENSIONAL PROBLEM (GRATING )USING 
THE SECANT-CORRECTOR SPECTRAL ITERATION APPROACH 
(MASTER'S THESIS BY ROBERT MIDDELVEEN) 



ABSTRACT 


The secant method is applied to an iterative algorithm 
of electromagnetic scattering from planar surfaces with 
periodic structure. The theory of convergent solutions for 
iterative techniques is discussed and examined. The Secant 
method is applied to the spectral iteration approach to 
accelerate and assure convergence of the basic iterative 
scheme. The derivation of the method as applied to surfaces 
containing parallel thin wire gratings. is presented, and the 
conditions for achieving convergence are explored. This new 
method is also applied to gratings made of coated wires. The 
reflection characteristics of the grating as a function of 
wire spacing, wire conductivity, and polarization of the 
incident field are computed, and the results are compared 
with those of previous works. Suggestions and recommen- 
dations for applying the method to more complicated struc- 
tures are also included. 
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I. INTRODUCTION 


Scattering from periodic structures such as grids or 
gratings has been of particular interest in the field of 
electromagnetics for many years. The field distribution 
about the structure caused by a plane incident wave, the 
induced current densities in the wires, and the reflection 
coefficient are the most important parameters used in 
designing grids and gratings. 

The objective of this thesis is to model the problem of 
electromagnetic wave scattering from a grating and compute 
the current density and reflection coefficient. Grids are 
important because they can be used as reflecting surfaces 
instead of solid metal surfaces, especially for modeling 
light-weight antennas for space applications. Moreover, by 
making use of the frequency dependence of these structures, 
they can be applied to filtering from the microwave to the 
optical wave regions. 

Many different methods have evolved for solving the 
problem of electromagnetic scattering from such structures. 
The most popular approach, the method of moments, usually 
requires large amounts of computer memory when applied to 
periodic structures. Another technique, the spectral- 
iteration technique (S.I.T.) developed by Tsao and Mitra [1] 
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circumvents this memory requirement/ but suffers from 
convergence problems. For example/ if the separation of 
adjacent wires/ or strips is less than two wavelengths/ then 
the S.I.T method will not converge. 

Brand [2] applied a corrective scheme that assured the 
convergence of the basic iterative equation for any wire 
spacing. This method, however, depends on the evaluation of 
numerical derivatives to generate a series of convergent 
iterations. In some cases the computation of the derivative 
can be so critical that the new corrective scheme fails to 
converge. This thesis presents an alternative, derivative- 
free technique which always converges for any spacing of 
adjacent wires, polarization of incident wave, and angle of 
incidence of the incoming wave. 

Another alternative method for solving scattering 
problems is the Fast Fourier transform-conjugate gradient 
method (FFT-C.G.) developed by Chistodoulou [3], This 
technique can be used to solve for either the strip currents 
or the electric field separately. Results obtained using 
this method are compared with those obtained using the new 
algorithm. 

Also included in this thesis is a study of electro- 
magnetic scattering from gratings made of coated wires. An 
internal impedance is used which takes into consideration 
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the effects of the substrate on the induced currents and 
reflection coefficient. This approach is particularly useful 
for space applications where a highly conductive coating is 
used in conjunction with a light substrate. Usually the 
electrical characteristics of the substrate are unimportant. 
For certain frequency ranges, however, the electrical fields 
and currents penetrate both coating and substrate so that 
the properties of both materials become important in the 
calculation of fields and currents. 



II. DERIVATION OF THE ORIGINAL ITERATIVE SCHEME 

A model is presented here that can determine the 
electric field and current density along the surface of a 


unit 

cell il 

lustrated in 

Figure 1. 
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magnetic potential 
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making use o 

f the free-space Greens func 

tion 

~. 

The v 

ector potential is 

a convolution of K 

and G, gi 

ven 

by 

F = 

/ G(r ' 

r’) K(F) dF' 




[2] 

where 

the f r 

ee-space Gre 

ens function is def 

ined by: 



G = 

(1/4 7T 

\~\ ) T Exp ( 

-j ) 



[3] 

The 

dyadic 

is denoted 

by T. The two vect 

ors F and 

k 

are 

i llus 

trated 

in Figure 

1. Returning to 

equation 

1, 

the 


magnetic field intensity as a function of magnetic vector 


4 



5 



the Plane 
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potential can be derived from Maxwell's equations and 
applying the Lorentz gauge condition ( see Appendix B ) . 

H* = -ju>£ ~ + (1/j Ufl ) V ( V • “) [4] 

Where /i is the permeability of the medium, and u is the 
angular frequency. Referring to Figure 1, there can be no 
magnetic current in the z direction because the planar 
structure is limited to the x-y plane. Therefore no 
component of the magnetic vector potential can exist in the 
z direction. Also because the stucture is located in the 
plane z = 0, G will be a function of x and y only. If the 
medium is allowed to be that of free space, then the 
permeability and permittivity will remain constant and are 
given by and £ 0 respectively. The propagation constant 

for this medium is defined by: 

k =(j [5] 

Equation 4 can now be expanded in the Cartesian coordinate 
system as follows: 

r s = l/j^J(kfF x + _^F X + d\ ) a x [6] 

L dx dv d x 2/ 


+ 

( kf Fy + e 2 F x + 

-£?yj a y] 


6x d y 

e i 2 J 


The subscript s in this equation signifies that this is the 
scattered field. The scattered field is caused by the 
incident field generating the magnetic currents in equation 
6. These currents in turn produce the scattered fields. The 
total H - field can be found by adding the incident and 
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scattered fields. In vector notation equation 6 is expressed 


= 1/j CJfic 


2 2 2 

K+d/dx 

d 2 /dxdy 


d/dxdy 
ko + d / dy 


Equation 2 is now substituted into equation 7. Taking the 
Fourier Transform of equation 7, the convolution of G and ~ 
in the space domain becomes a multiplication in the Fourier 
domain. The transformed scattered magnetic intensity is 


given by: 


H s = l/jo>// c 


,2 2 
k Q ~ a mn 




a mn^mn 

.2 a 2 
k 0 “ P mn 


G K 


The tilde symbol is used to denote the transform of the 
variable in the Fourier domain. The parameters a mn ' and 
f} mn are referred to as the Floquet modes and they are 
defined as follows, reference [2]: 


Ojnn =27 T m/a 


- k 0 sin $ cos 0 


Amn = 27rn/c - (27rm/a) cot£ - k 0 sin# sin0 


[ 10 ] 


Their values depend upon the cell geometry of the planar 
surface being studied. The angles $, <p , and Q are shown 
in Figure 1. The number of sampling points across the unit 
cell in the x and y direction are given by m and n 
respectively. The Floquet modes allow for the effects of 
coupling between the conducting regions of the planar 
surface. The Fourier Transform of Green's function is given 


G - - ( j / 2 ) ( k„^ - O mri — fimn) 1 
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If the Inverse Fourier Transform is applied again to 
equation 8, the scattered field in the space domain will be 


given by: 


k° - ®mn tt mn Ann ^ 

H s = l/j(jp 0 J2 G K Ex P(j[«mn x [12] 

mn -ot mnAmn Amn 

By using the equivalence theorem and applying the 
appropriate boundary conditions to the scattered H field at 
z - 0 r the total tangential H field .can be solved in terms 
of the transformed electric field in the aperture as: 

_ F ®mn Ann ®mn 1 S — 


H tinc = G E Exp(j[a mn x +£ mn y]) [13] 
mn -lc o + Amn -tt mn Amn 

To include the contribution of the "h field along the 
conducting strip, the current densities have t-o be added to 
equation 12 to yield, reference [2]: 


®mn Amn ** mn £■; ~ 

Trc'(J) = H tinc + 2/jOJUoZ G E C 14 ] 

mn ~^o + Amn ~®mn Amn 
Ex P ( j [®nm x + Amn Y ] ) 

Because the current density can only be present on the 
conducting strips, the truncation operator Trc and its 
complement are introduced. These are defined by: 


Ter [ Xfr) ] = 


X(r) for r in the aperture 

I 

0 for r in the conducting region 


[15] 


Ter ' [ X ( r ) ] = 


0 for r in the aperture 

1 

X(r) for T in the conducting region 


[16] 
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In equation 14, direct solution for the electric field is 
not possible since both the strip current and electric field 
are unknown. For this reason Tsao and Mittra [4] developed 
an iterative equation to solve for both the electric field 
and the strip current. Returning to equation 14, the 
following simplification is made: 

~ f Ofon/Smn C*nin 1 ~ 


" k o + Ain ®mn^mn 
With this substitution, and th'e fact that the tangential 
field is present only in the aperture and the current 
density exists only along the conductor, equation 14 can be 
written as: 


[17] 




Trc’(J) = T rc'(ir tinc + 2/jUtflo F 1 [G 2 F ( Trc [ E t ] ) ] ) [18] 

Similarly the tangential electric field can be derived from 
equation 14: 

“ t = F _1 [G 2 -1 F ( / 2 [Trc * [7] - *H t inc l) ] [19] 

This electric field represents the field across the entire 
cell. The field is also valid on the conducting strip 
because the truncation operation is performed on this field 
in equation 18. Equation 18 is substituted into equation 19 
yielding the basic form of the iterative equation in terms 
of the electric field. 

F ti = F _1 [g" 2 _1 F ( j / 2 [Trc' (ir tinc + 2/ju>v o F _1 

[G 2 F (Trc [E^.! ])]) - H tinc ])] 


[ 20 ] 
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Brand [2] imposed his corrective technique on this 
equation. As an alternative to Brand's correction, the 
secant technique can also be applied to equation '20. This 
technique will eliminate the convergence problems caused by 
the unavailability of an analytic derivative. 



III. THE SOLUTION OF FIXED POINT PROBLEMS VIA 
ITERATION FUNCTIONS 

General Theory of Iterative Functions 
Iterative functions are often constructed to solve 
equations of the following form: 

f ( x ) = 0 [21] 

and F(“) = 0 [22] 

Equation 21 represents a single variable complex function. 
In equation 22, both the function and argument are vector 
quantities whose elements can be complex. If these 

functions are of a very complex nature, a direct solution 
may not be available. However, a solution of arbitrary 
precision may be obtained using iterative techniques. These 
techniques make use of equations 21 or 22 to construct 
iterative functions of the form: 

x i+l = 9( x i ) [23] 

Xi +1 = GUi) [24] 

Equation 23 is a single variable function and equation 24 is 
the vector equivalent. 

Solutions of f(x)=0 

The iterative process consists of starting with an 
initial guess x 0 , and inserting this value in equation 23 to 


11 
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obtain a new iterate x^. The process is repeated until 
x i+1 = x^ + e where e is an allowable error. If x*= x i+1 = x i , 
then x* is called the fixed point or the best numerical 
solution of equation 23. The fixed point of g(x) is a root 
of f(x). There are two conditions required to assure a 
convergent solution: 

1) On a closed region containing the solution x, g(x) 
should be continuous. 

2) For any arbitrary points s and t in this region the 
following condition must be met: 

| g (s ) - g (t )| < p | s — 1 1 [25] 

0 < p< 1 

These conditions imply that g(x) must be differentiable 
over the interval of interest and that the magnitude of its 
derivative must be less than unity. Froberg [4] has an in- 
depth proof of this statement. If the above conditions hold 
then g(x) is said to be a contraction, and the iterative 
process will eventually produce a fixed point x*. The 
graphical description of the iteration process is 

illustrated in Figure 2. The fixed point x* is obtained by 
finding the intersection of y = x and y = g(x). The first 
four iterations are included beginning with the initial 
guess denoted x c . 
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Figure 2. Graphical Illustration of the Iterative Process. 
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To illustrate the iterative process, consider the 
following equation: 

0 = f ( x ) = + 1 1 x 3 + 35. 5x 2 + 57x 1 + 31.5 [26] 

This equation contains roots at -1, -7, -1.5-jl.5 and 

-1.5+jl.5. It is desired to find these roots using an 
iterative process. An obvious choice for g(x) is obtained by 
solving equation 26 for x 1 . 

g ( x ) = (-l/57)(x U + llx 3 + 35. 5x 2 + 31.5) [27] 


The derivative of this function is: 
g ' ( x ) = ( -1/57 ) ( 4x 3 +33x 2 +71x) 


[28] 


The values of this derivative at the different roots are: 


g ' (-1) = 0.7368 [29] 
g < (- 7 ) = 4.4211 [30] 
g ' ( —1 . 5— j 1.5) = 1.4193 Exp ( -j 0 . 1865 ) [31] 
g ' ( -1 . 5+ j 1 . 5 ) = 1.4193 Exp(j0.1865) [32] 


Note that f(x) and g(x) are both continuous over the entire 
complex plane but that the derivative of g(x) is less than 
unity only in the interval about the root at x = -1. With an 
appropriate initial guess g(x) should therefore converge to 
the fixed point x* = -1. The first 20 iterations, generated 
by equation 27, when x„ = 2 are tabulated in Table 1. It 
should be mentioned that this choice of g(x) cannot be used 
to find the other roots, and for an inappropriate initial 
guess the iterative process may diverge altogether. 
Fortunately, there are other techniques which can be applied 



15 


TABLE 1 

EXAMPLE OF ITERATIVE EQUATION DERIVED FROM f(x). 
CONVERGENT ABOUT THE REGION x = -1.0 


Iteration 

9 ( X-) 

0 

-4.868421 

1 

-2.901527 

2 

-2.325328 

3 

-2.006732 

4 

-1.785657 

5 

-1.618081 

6 

-1.485960 

7 

-1.380177 

8 

-1.295302 

9 

-1.227566 

10 

-1.174003 

11 

-1.132097 

12 

-1.099659 

13 

-1.074794 

14 

-1.055895 

15 

-1.041630 

16 

-1.030925 

17 

-1.022926 

18 

-1.016970 

19 

-1.012546 


f (x) 
g (x ) 


x U + llx 3 

^1 { x^ + 
57 


+ 35. 5x^ + 57x 
llx 3 + 35. 5x^ 


+ 31.5 


+ 31.5 


) 
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to this problem to assure convergence and increase the 
convergence rate. 


Techniques to Accelerate and Assure Convergence 

The techniques used to resolve these problems involve 
the synthesis of a better iterative function than the 
example illustrated in equation 27. These methods involve 
the use of derivatives of the function f(x) and the 
knowledge of previous iterates stored in memory. 


One of the techniques that can be used to find 
roots of equation 26 is the Newton-Raphsom method, 
iterative function synthesized with this technique 
the use of derivatives of f(x) and is written 
following form: 

g(x i ) = - f (x t )/f ' (x± ) 


all the 
The new 
involves 
in the 


[33] 


This new iterative function has a derivative with a 
magnitude that is always less than unity over the complex 
plane if f(x) and f'(x) are well-behaved functions. The 
problem encountered with the previous sample iteration 
function is no longer present and equation 33 can now be 
used to solve for the missing roots. Table 2 lists the first 
10 iterations of the solutions to the real roots of equation 
26. From this table it can be seen that the rate of 
convergence has also been increased. This acceleration is to 
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TABLE 2 

EXAMPLE OF THE NEWTON-RAPHSON I.F. APPLIED TO f(x) TO 
ACCELERATE CONVERGENCE AND FIND ADDITIONAL ROOTS 


Iteration (i) 

x 0 = 2.0 

g (>4 ) 

x 0 = -10.0 

g Ui ) 

0 

0.921488 

-8.513304 

1 

0.117348 

-7.575745 

2 

-0.470798 

-7.121042 

3 

-0.846059 

-7.006845 

4 

-0.985979 

-7.000025 

5 

-0.999888 

-7.000000 

6 

- 1.000000 

-7.000003 

7 

- 1.000000 

-7.000003 

8 

- 1.000000 

-7.000002 

9 

- 1.000000 

-7.000000 


f ( x ) = x k + 1 1 x 3 + 35. 5x 2 + 57x + 31.5 


= x i - 


g ) 


f (Xi )/f ' (Xi ) 
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be expected as more information about f(x) is used in the 
formulation of new iterations. The technique is also valid 
for complex valued functions. The first 10 iterations 
solving for the complex root are shown in Table 3. A 
graphical illustration of this technique is depicted in 
Figure 3. The first three iterations are shown beginning 
with the initial guess x 0 . The fixed point is denoted by 
x* . Note that the true tangent is used to calculate the 
next iterate. 

The Newton-Raphson method is universally known and is 
the most popular and useful iterative function. 
Occasionally, when an analytic derivative is not available, 
one can be approximated by perturbing the original function 
by a small increment, denoted by del ( 4 ). However, the error 
introduced using this method may be critical. If del is too 
large, the approximation will not be valid at the desired 
point. Figure 4 shows how the numeric derivative will 
contain a large error if the function is changing too 
rapidly with x. If del is too small the approximation is 
limited by the precision of the numerical operation. 

The secant method does not depend on the evaluation of 
any numerical derivative. The secant iterative function is 
defined by: 

g(x i ,x i _ 1 ) = x i - f(x i )[(x i -x i _ 1 )/(f(x i )-f(x i _ 1 )] [34] 
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TABLE 3 

EXAMPLE OF COMPLEX VALUED NEWTON-RAPHSON I.F. APPLIED TO 
f ( x ) TO ACCELERATE CONVERGENCE AND FIND COMPLEX .ROOTS 


Iteration (i) 

Xo = 

-4.000000 - 
g(x i ) real 

j4. 000000 
g ( x^ ) imaginary 

0 


-3.485554 

-2.714932 

1 


-2.896810 

-1.753410 

2 


-2.228628 

-1.196201 

3 


-1 . 583420 

-1.059435 

4 


-1.152799 

-1.751122 

5 


-1.370047 

-1.506637 

6 


-1.495772 

-1.483877 

7 


-1.500129 

-1.500248 

8 


-1.500000 

-1.500001 

9 


-1.499999 

-1.500000 


f (x) 


x^+ 1 1 x 3 + 35. 5x 2 + 57x + 31.5 


g (x i ) = x 


f ( )/f ' ( x ± ) 
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an 


Estimated tangent 


21 



rivative 


22 


This method is based on a similar geometry as that of the 
Newton-Raphson method. Figure 5 shows that sequential 
iterates are used to estimate the tangent. This method tends 
to converge more slowly but always locks in to the fixed 
point. The previous method/ when the derivative was 
approximated numerically/ occasionally oscillated about a 
point in the vicinity of the solution. The first 10 
iterations of the solutions for the real roots of the sample 
function are shown in Table 4. The first 10 iterations 
solving for the complex root are shown in Table 5. As with 
the Newton-Raphson method/ the secant method converges much 
more rapidly than the original iterative function given by 
equation 27 because more information is being used. In 
addition, convergence is assured for all roots. 


Formulation of the Problem 

The preceeding example was concerned with single valued 
complex functions, such as those described by equation 1. In 
some useful applications, however, the domain and range of 
the function are vector quantities such as those described 
by equation 22. Such is the case with the original iterative 
function derived by Tsao and Mittra [1]. This function (see 
chapter II) is repeated below for convenience: 

T ti = F -1 [ G~F ( j(Jfl 0 /2 [ Trc ' ( H tin J 2/jU)U 0 F _1 
[ G 2 F (TrctE”^]) - H tinc ])] 


[ 20 ] 
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TABLE 4 

EXAMPLE OF THE SECANT I.F. APPLIED TO f(x) TO ACCELERATE 
CONVERGENCE AND FIND ADDITIONAL ROOTS 



x_i = 3.0 

x — i = -11.0 


Xg =2.0 

o 

• 

0 
1 — 1 

1 

II 

o 

X 

Iteration (i) 

g (*i ) 

g (*i ) 

0 

1.230089 

-8.808706 

1 

0.587572 

-8.044291 

2 

0.068331 

-7.486720 

3 

-0.351034 

-7.169444 

4 

-0.671497 

-7.034141 

5 

-0.881845 

-7.002733 

6 

-0.977097 

-7.000047 

7 

-0.998424 

-7.000005 

8 

-0.999980 

-6.999999 

9 

-1.000000 

-7.000000 

f ( X ) = x 1 * + llx 3 

+ 35. 5x^ + 57x 

+ 31.5 

9< x i' x i-l > - x i - 

f ( x. ) [ ( x.-x. . 

i l i-l 

)/ ( f ( x.^ )-f (x i _^ ] 
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TABLE 5 

EXAMPLE OF COMPLEX VALUED SECANT I.F. APPLIED TO 
f ( x ) TO ACCELERATE CONVERGENCE AND FIND COMPLEX ROOTS 



x_i = 5.000000 + 

j5. 000000 


Xq sz 4.000000 + 

j4. 000000 

Iteration (i) 

g ( x . ) real 

l 

g(x ) imaginary 
i 

0 

-3.827079 

-3.051936 

1 

-3.365816 

-2.280450 

2 

-2.940393 

-1.648176 

3 

-2.441053 

-1.230939 

4 

-1.956624 

-1.021398 

5 

-1.482012 

-1.010186 

6 

-0.944583 

-1.453021 

7 

-1.845381 

-1.300908 

8 

-1.660885 

-1.734051 

9 

-1.567421 

-1.431206 

f ( x ) = x U 

+ llx 3 + 35. 5x 2 + 57x 

+ 31.5 

g (*i / x i _ 1 ) 

= x.^ - f ( x ± ) [ ( x i - x i _ 

)/( f ( x i ) - f ( x i _ 1 ) ] 
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To simplify the above expression/ the left side of the 
equation 20 is defined as an operator of the electric field 
vector. 


E. = L ( E. ) 
1 1 l-l 


[35] 


The iterative equation described by equation 20 does not 
always converge to a solution for the electric field. This 
is analogous to the convergence problems encountered with 
the sample iterative function, equation 27. Specifically, 
equation 20 was found only to converge for very large wire 
spacings. There are two equivalent techniques to alleviate 
this problem. Brand [2] chose to relax equation 20 using the 
following equation: 

e x - v r i-i> - r e i-i + < 1 - r > Wi > t 36 ’ 

where 1„ and e are the individual elements of the l”_ and e" 
2 2 

vectors resectivly and 


R = 1 1 ' (E) / ( 1 1 * (E) - 1 ) [37] 
It is in equation 36 that problems with convergence emerge. 
This equation requires the derivative of . Because an 
analytic derivative is not available, an approximate 
derivative is defined by: 


1 '( E ) = lil+4 ) - KB) [38] 

A 

It was this approximation that was often found to be 
inadequate. Because of intrinsic differences between the IBM 
and VAX computers and their compilers, Brand's results could 
not be duplicated on the VAX due to the limitations in 
precision when approximating the derivative. 
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An alternative approach was taken here/ by defining a 
new equivalent vector function, which has as its root the 
solution of the electric field. This new function is defined 
by: 

F (If) = L 1 (E) - E [39] 
The Newton-Raphson and Secant techniques can be applied 
directly to equation 39. The vector F is called the residue 


vector and has a value proportional to the remaining error 
in the electric field. Table 6 compares the value of this 
vector after 10 iterations for the S.I.T. with the 
contraction and secant correctors applied. It can be seen 


that the secant method produces a much smaller residue. The 
Newton-Raphson method can be applied to equation 39 as 
shown below: 

1 2 (E) = e - f ( e") / f ' (E) [40] 
where f represents an individual element of E 7 . Equations 36 
and 40 are mathematically identical. Brand [2] included a 
formal proof repeated in Appendix C showing this 
equivalence. The Newton-Raphson iterative function produced 
results virtually identical to Brand's original model. 
Convergence problems persisted for certain input parameters 
such as low angles of incidence and small wire spacings. 

The secant method can be applied to equation 14 in the 


following fashion: 
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TABLE 6 

COMPARISON OF RESIDUE VECTORS 


Array 

Element 

S.I.T. 

(Real) (Imaginary) 

S.C.S 
( Real ) 

. 1 . 

( Imaginary 

1 

0.003697 

0.004007 

-.00000054 

0.00000182 

2 

0.003874 

0.002393 

-.00000024 

-.00000909 

3 

0.003956 

0.001639 

0.00000703 

-.00001219 

4 

0.004024 

0.001026 

-.00000775 

-.00000885 

5 

0.004069 

0.000608 

-.00000519 

-.00000742 

6 

0.004109 

0.000244 

-.00000465 

-.00000682 

7 

0.004139 

-.000029 

-.00000477 

-.00000659 

8 

0.004166 

-.000273 

-.00000471 

-.00000641 

9 

0.004168 

-.000459 

-.00000477 

-.00000635 

10 

0.004204 

-.000625 

-.00000477 

-.00000632 

11 

0.004218 

-.000749 

-.00000477 

-.00000626 

12 

0.004230 

-.000857 

-.00000483 

-.00000626 

13 

0.004238 

-.000930 

-.00000471 

-.00000626 

14 

0.004244 

-.000989 

-.00000477 

-.00000620 

15 

0.004247 

-.001018 

-.00000477 

-.00000626 

16 

0.004249 

-.001033 

-.00000471 

-.00000626 

17 

0.004247 

-.001018 

-.00000471 

-.00000626 

18 

0.004244 

-.000989 

-.00000471 

-.00000626 

19 

0.004238 

-.000930 

-.00000483 

-.00000620 

20 

0.004230 

-.000857 

-.00000477 

-.00000608 

21 

0.004218 

-.000749 

-.00000477 

-.00000620 

22 

0.004204 

-.000625 

-.00000471 

-.00000620 

23 

0.004186 

-.000459 

-.00000471 

-.00000629 

24 

0.004166 

-.000273 

-.00000471 

-.00000641 

25 

0.004139 

-.000029 

-.00000471 

-.00000656 

26 

0.004109 

0.000244 

-.00000477 

-.00000679 

27 

0.004069 

0.000608 

-.00000519 

-.00000739 

28 

0.004024 

0.001026 

-.00000781 

-.00000888 

29 

0.003956 

0.001639 

0.00000703 

-.00001213 

30 

0.003874 

0.002393 

-.00000024 

-.00000909 

31 

0.003697 

0.004007 

-.00000060 

0.00000188 

32 

0.003328 

0.007377 

0.00000649 

-.00001098 
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1,(E . ,E. . ) = e . - f ( E . ) ( e , - e, ) 

3 1 l-l l li l-l 

where 1^ is an individual element of 


operator L ^ . 


/ f(E*. )-f( 

l 

the vector 


E. n ) [41] 

l-l 

produced by 


Although the vector and single valued iterative 
functions appear to be constructed in identical fashion, 
there are important differences regarding the number of 
solutions and the conditions necessary to assure 
convergence. The vector operator L is called a contraction 
operator in a particular domain if it satisfies the 
following condition: 

d 2 [L(U) ,L(V) ] p d 2 [U,V] [42] 
The magnitude of p is always less than unity and TT and T are 
any two vectors in the domain. The operator d 2 referred to in 

equation 42 is the distance function and is defined by: 

__ _ fn o 1 1/2 

d 2 (U ' V) = | l u i" v il j [43] 

If L is a contraction throughout a given domain, then from 
any starting point within that domain, there will be one and 
only one fixed point defined by: 

TT*= L(TT*) [44] 
The formal proof and a discussion of vector spaces is 
provided by Stakgold [5]. 


Brand [2] proves with mathematical rigor that the 
Newton-Raphson iterative equation, given infinite numerical 
precision, will always converge to the fixed point. 
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Restrictions are imposed that are even more stringent than 
those of equation 42. It is shown that the Newton-Raphson 
method complies with these new conditions. Given that the 
Secant method is an approximation of the Newton-Raphson 
method, it is intuitive that this method should also force 
convergence upon the original iterative scheme devised by 
Tsao and Mitra. In fact Traub [6] shows that the order of 
the iterative function for the secant method is 1.62 as 
compared with 2 for the Newton-Raphson method. Any two 
vectors in the domain of L can be chosen as the U and V of 
equation 42. Brand [2] chose as his two vectors E and 
E - + A . To monitor the performance of the contraction of 
the secant iterative function, Ej^ and E i+1 proved to be 
convenient vectors. A contraction factor is defined below 
for the secant method: 

Con = d(L(Ej ),L(j i _ 1 ) [45] 

dT^TST^) 

This factor was verified to comply with equation 17 when the 
contraction process was taking place. Applying the secant 
method, it was found that convergence was obtained for any 
wire spacings, any polarization of the incident wave, and a 
wide range of wire conductivity. 



IV. THE INTERNAL IMPEDANCE OF THE WIRES IN THE GRATING 


To account for the finite conductivi 

the following boundary condition must be 

E . . + E = Z I 

tine s 

where Z is the internal impedence of the 
the current present in the conductor. It 
derive an expression for the impeadance. 


ty in coated wires, 
met : 

[46] 

conductor and I is 
is now necessary to 


At very high frequencies the impedance of a solid wire 
can be obtained using the impedance formulas for a semi- 
infinite plane solid. At a sufficiently high frequency the 
curvature of the wire becomes unimportant. This occurs when 
the skin depth for the conductor becomes small compared with 
the radius. The wire may then be considered a plane solid 
with infinite depth and a width equal to its circumference. 
The internal impedance of the wire for this case is given 
by : 

Z = Z / 2 7T r [47] 
s 

where Z is the impedance of the plane solid, and 2 n r is 
the circumference of the wire. Z is then measured in Ohms 
per unit length. Z g for a good conductor is given by 
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Z s = ( 1 + j ) = R s ( 1 + j) [48] 

Ob 

where a is the conductivity of the conductor in Siemens/ A 
is the skin depth in meters and Rs is the surface 
resistivity in ohms per square. The skin depth, 6 , is: 

d = 1 [49] 

y[ muo 

where f is the frequency, and n is the permeability of the 
conductor. The surface resistivity therefore becomes: 

R s = = JnfM/o [50] 

ob y 

Although the actual mesh structures of interest are 

made of small round wires, the algorithm used to estimate 

the characteristics of this mesh apply only to a planar 

structure of negligible thickness. Therefore the conductors 

have to be modeled as conducting strips. For scattering 

problems this modeling is done using the concept of 

equivalent radius. The wire of radius r can be replaced by a 

strip of width w, where r is given by: 

r = 0.25w [51] 

The final expression for the impedance of the wire using 

equations 47, 48, and 51 is given by: 

Z = 2R S _U + j) [52] 

71 w 

Finally the mesh is usually not made of solid wire but a 
coated material. In this way full advantage can be taken of 
both the mechanical characteristics of the substrate and the 
electrical characteristics of the coating. If currents and 
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fields are able to penetrate both materials, then the wire 
impedance will no longer be predicted by equation 51. A more 
accurate measure of impedance, [7] derived in Appendix C, is 
given by: 

Z = 2 R g (1 + j)[ sinh(T^d) + (Rg^/Rpi )cosh(T 1 d) ] [53] 
n w [ cosh(T^d) + tRg2 /^l ) s i nh ( Tf d ) ] 

where, = ( 1 + j ) = (1 + j ) ^ 7T f [54] 

d is the thickness of the coating in meters, and Rg^ and R^ 
are the surface resistivities of the surface and the sub- 
strate respectively. For very small values of coating 
thickness, equation 53 reduces to the impedance of a wire 
made of only the substrate. For large values of thickness 
the wire appears as if it is made entirely of the coated 
material. For intermediate values of thickness the resis- 
tive and reactive parts of the impedance are no longer equal 
in magnitude. Figures 6 and 7 show the wire impedance, 
normalized with respect to that of the coating material, for 
varying ratios of thickness to the skin depth of the 
coating. Figure 6 corresponds to a solder coating on a 
copper substrate, with %2 /%l = 0.34. Figure 7 corresponds 
to a silver coating on a brass substrate, with 


%2 " 1 * 6 * 






V. RESULTS 


First the fields and currents, as calculated by the 
secant method, will be examined and compared with other 
published results. There are two sources of published 
results with which the new algorithim is compared. Brand [2] 
used the spectral domain approach with the contraction 
factor denoted by S.I.T, and Chistodoulou [3] used the the 
fast Fourier transform-conjugate gradient method denoted by 
FFT-C.G. 

An important parameter that needs to be mentioned at 
this point is the number of sampling points required to 
represent the physical situation. For very thin wires a 
greater number of sampling points will produce more accurate 
results because there will be less quantitization error of 
the strip width. The position in the cell at which each 
sample is taken can effect the results. In addition, the 
greater the number of sampling points the slower the 
contraction process will be. 

Table 7 shows a comparison of the three methods for a 
grating of very thin wires. For each case the current 
density of the wire was computed. All cases were examined 
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with the incident field normal to the plane and the E field 
parallel to the wires. There were 32 sampling points for 
each unit cell with one point laying on the strip. The 
results of these methods are in good agreement. Any of these 
methods can be used to predict the current density on a 
strip for different wire spacing and wire thickness. 

Figure 8 illustrates the current density across a wide 
strip. Sixteen sampling points lie on both the strip and 
the aperature. The incident field generating this current is 
again normal and copolar. The current density is seen to be 
very large at the edges of the strip. This result shows that 
the new algorithim can predict edge effects. 

Table 8 shows a comparison of the electric fields 

across the entire cell as predicted by the three methods. It 

should be mentioned that the FFT-C.G. method does not 

actually compute the field across the strip region, it is 

assumed zero because the conductivity of the strips is very 

large. There are again 32 sampling points with two points 

lying on the strip. Note that the electric field located on 

the strip is predicted to be much lower for the S.C.S.I. 

method than for the S.I.T. For perfectly conducting strips, 

low field values on the strip indicate that the boundary 

condition IT . + E = 0 is satisfied. This condition 

tine s 

could always be used as an accuracy-check. 
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TABLE 8 

COMPARISON OF THE ELECTRIC FIELD IN A UNIT CELL 


CELL POINT 


-0.129126310 

-0.120795548 

-0.112464786 

-0.104134083 

-0.095803320 

-0.087472617 

-0.079141915 

-0.070811152 

-0.062480479 

-0.054149747 

-0.045819018 

-0.037488285 

-0.029157557 

-0.020826824 

-0.012496091 

-0.004165362 

0.004165362 

0.012496091 

0.020826824 

0.029157557 

0.037488285 

0.045819018 

0.054149747 

0.062480479 

0.070811152 

0.079141915 

0.087472618 

0.095803320 

0.104174083 

0.112464786 

0.120795548 

0.129126310 


S.I.T. 


0.182258561E-1 

0.384013295 

0.560089946 

0.661797166 

0.738476992 

0.796682596 

0.844159245 

0.882627010 

0.914659142 

0.940926552 

0.962553382 

0.979851842 

0.993371725 

1.00326443 

1.00979042 

1.01301479 

1.01301575 

1.00979042 

1.00326347 

0.993371725 

0.979854842 

0.962553859 

0.940926552 

0.914659023 

0.882627010 

0.844159365 

0.796682477 

0.738476753 

0.661787643 

0.560090780 

0.384014010 

0 . 1 82 278380E-1 


FFT-C.G. 


0.000000000 
0.377430677 
0.556849957 
0.660458642 
0.738587141 
0.797879934 
0.846264482 
0.885471702 
0.918128848 
0.944895148 
0.966925740 
0.984547973 
0.998328090 
1.00839138 
1 .01503468 
1.01830196 
1.01831055 
1.01503181 
1.00839233 
0.998323321 
0.984559417 
0.966914829 
0.944902539 
0.918127894 
0.885474324 
0.846258521 
0.797884822 
0.738585949 
0.660460711 
0.556844115 
0.377436161 
0.000000000 


S.C.S.I. 


4 . 8894606E-5 

0.3747891 

0.5529166 

0.6556744 

0.7332161 

0.7919211 

0.8399065 

0.8787327 

0.9110555 

0.9375529 

0.9593651 

0.9768083 

0.9904386 

1.000410 

1.006989 

1.010238 

1.010238 

1.006989 

1.000410 

0.9904386 

0.9768084 

0.9593652 

0.9375528 

0.9110555 

0.8787327 

0.8399065 

0.7919205 

0.7332161 

0.6556743 

0.5529164 

0.3747890 

4. 8571634E-5 
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Figure 9 shows the electric field across a unit cell 
with a large strip size. The electric field magnitude along 
the strip is seen to be very small. The shape of the 
electric field is in good agreement with results published 
by Brand [2], 


Next/ the reflection coefficent predicted using the 
secant method is compared with FFT-C.G and S.l.T. methods. 
In addition/ comparisons are made with results published by 
Wait [8]. The reflection coefficient for cell widths less 
than one-half wavelength is equal to the first element in 
the transformed electric field vector, i.e., the first mode. 


r = e~ = eT„- e. . 

1 00 00 tine 


[55] 


For these cell widths only one propogating mode will be 
present so that only one array element is needed. These 
narrow cell widths are important because when the wire 
spacings become 0.5 wavelengths or smaller, the planar 
surface begins to resemble a solid reflector. 


In Table 9 the magnitude of the reflection coefficient 
for different values of wire spacing for normal incidence 
and wire radius of 1/600 wavelengths is compared with other 
methods. As expected, the grating begins to appear as a 
solid reflecting surface as the wire spacing becomes small. 
All three methods are in very good agreement. 
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TABLE 9 

COMPARISON OF REFLECTION COEFFICIENTS 
FOR DIFFERENT WIRE SPACINGS 


SPACING FFT-C.G. S.C.S.I. S.I.T 


0.125 

0.844 

0.10 

0.888 

0.06 

0.954 

0.05 

0.967 

0.02 

0.994 

0.01 

0.999 


0.844 

0.843 

0.892 

0.885 

0.971 

0.960 

0.985 

0.969 

0.999 

0.994 

1.000 

0.999 




44 


Figure 10 depicts the changes in the reflection coef- 
ficient for cell widths of 1/2, 1/4, and 1/8 wavelengths as 
0 is varied from 0 to 90 degrees. The wire radius remains 
constant at 1/600 wavelength. The angle <p remains constant 
at zero degrees corresponding to TE polarization of the 
electric field. These results are in good agreement with 
results published by Wait [8} and Brand [2], 

Figure 11 illustrates the behavior of the reflection 
coefficient when the angle <t> is held constant at 90 
degrees. The cell width is held constant at 1/4 wavelength 
while $ is varied from 0 to 90 degrees. The wire radius is 
again 1/600 wavelength. For this angle of <f> it is observed 
that there exists a region of maximum transmission at 9 
equal to approximately 67 degrees. This is analogous to the 
Brewster angle associated with dialectric materials. 

The effects of the substrate conductivity of coated 
wires on the reflection coefficient were studied next. In 
Figure 12 the conductivity of the coating was held constant 

O 

at 5(10°) Siemens while the substrate remained constant at 
50 Siemens. The top curve corresponds to a very thick 
coating so that the conductivity of the strips is equal to 
the conductivity of the coating alone. The lower curve 
corresponds to a very thin coating so that the conductivity 
of the strip is greatly reduced. The complete set of curves 
illustrate the effects of varying the coating thickness. It 




Figure 10. Reflection Coefficient 
for Varying Values of Cell Width and(J! 
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Angle . 
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Figure 12. Effects of Finite Surface Conductivity and Coating Thickness. 
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can be seen that as the thickness of the coating becomes 
smaller, the incident field penetrates deeper into the sub- 
strate. Since the substrate is of lower conductivity than 
the coating, it is expected to have more current losses in 
the wire and hence a lower reflection coefficient is 
obtained . 

Figures 13, 14, 15 and 16 depict the effects of the 

substrate alone. The conductivity of the coating is held 

constant at 10 Siemens while the substrate is allowed to 

g 

vary from 10 to 10 Siemens. Note that the x axis cor- 
responding to the substrate conductivity has a log scale for 
each case. When the angle of incidence is perpendicular to 
the plane, the substrate is seen to have the most pronounced 
effect as shown in figures 13 and 14. The reflection coef- 
ficient increases as expected with an increase in strip 
conductivity. As the angle of the incident wave becomes much 
lower in" relation to the plane, the reflection coefficient 
remains more constant as the substrate conductivity changes 
as shown in figures 15 and 16. An interesting phenomenon 
appeared when <p and $ were set equal to 90 and 70 res- 
pectively. These particular values of 0 and $ are 

analogous to the Brewster angle of dialectrics. At this 
angle, the grating no longer behaves entirely as predicted 
for a reflecting surface. As the overall strip conductivity 
increases, the reflection coefficient drops slightly. At 
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this angle the surface appears more as a dialectric material 
with its permeability changing. 

These results were all run on the VAX 11/750. For 32 
sampling points with 0 equal to 0 degees and Q varying, 
the secant method converged to the third decimal of the 
reflection coefficient within 8 to 12 iterations. For very 
low angles, 6 being greater than 88, more iterations are 
needed. Convergence took 290 iterations when Q was equal 
89. There is also a correlation between the number of 
iterations and the number of sampling points. Table 10 shows 
that convergence is more difficult to obtain with a greater 
number of sampling points. 


TABLE 10 

SAMPLING POINTS VS. ITERATIONS REQUIRED 


Sampling points 


Iterations 


32 

64 

128 

256 

512 


8 

8 

16 

14 

20 


The 32 sample points with 8 iterations require 5.25 
seconds of CPU time while 512 sampling points with 20 
iterations requires 38.58 seconds. 





VI. SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS 

Brand's contractor corrector method failed to converge 
on occasion for various input parameters when running on the 
VAX computer. The reason for this failure was found to be 
the error introduced by using a numeric estimate of the 
derivative . 

An alternative derivative-free method was developed to 
insure the convergence of the spectral iteration approach as 
applied to the electromagnetic scattering from gratings. 
This method was derived by beginning first with the original 
spectral-iterative equation to which no correction was made. 
The general theory of iterative techniques was then covered. 
Basic examples were presented illustrating how the secant 
technique could be applied to solve single-valued complex 
functions. Finally, the secant method was applied to the 
vector space used by the spectral-iteration equation. 
Alternative methods for solving electromagnetic scattering 
are always of interest because the criteria for obtaining 
solutions become more stringent as the geometry of the 
problem becomes more complex. 

This new method was used to solve for the currents and 
fields lying in the plane of the grating. The reflection 
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coefficient was computed for cell widths under 1/2 
wavelength. These results were compared with previously 
published data and found to be in good agreement. 

The solution of scattering from a grating is a one- 
dimensional application of the spectral iterative technique. 
The cell geometry only changes along one axis of the plane 
containing the grating. A grid would require a two- 
dimensional application of the basic technique because the 
geometry changes along both the x and y axis. The contractor 
corrector has been applied to this two dimensional 
configuration and failed to yield a convergent scheme. 
Because the conditions necessary to assure convergence are 
more stringent, the error introduced by using a numeric 
derivative could be the critical factor causing this 
failure. The secant method has yet to be applied to this 
type of problem and could possibly provide a method of 
solution leading to convergence. 

No spectral iterative method has as yet been applied to 
geometries more complex than grids or gratings. For the one- 
dimensional problem the cell could contain various strips of 
varying size. The two-dimensional problem can be that of 
virtually any repeating planar structure. This case is 
important because surfaces approximating reflectors are not 
usually a grid, but a mesh structure which can have a very 
complex geometry. The spectral iterative techniques are 
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particularly well suited for these types of studies 
only the truncation operators are geometry dependent 
recommended that wire mesh geometries be studied to 
that convergence will still take place and that the 
remain acceptable. 


because 
It is 
verify 
results 



APPENDIX A 


MAIN PROGRAM AND ASSOCIATED SUBROUTINES 

The following program is written in FORTRAN 77 source 
code. It should run on most FORTRAN compilers and has run 
successfully on the VAX and IBM PC computers. This program 
will solve for the electric field across the aperature of a 
unit cell consisting of parallel wires. This unit cell is 
the repeating section of an infinite grating. The program 
will also solve for the current densities on the wires and 
the reflection coefficient of the grating. The program is 
presented in its interactive version, with appropriate 
prompts to request input data. This appendix consists of a 
summary of variables, subroutines and a program listing. 
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HI Incident magnetic field 

Z Internal impedance of the wires 

CREF Reflection coefficient 


CONVERGED Boolean indicator of convergence 
K Propagation constant 

CK_ Constants used in the original iterative eq. 

RS1/RS2 Surface resistivities of coating and substrate 
RATIO Ratio of coating thickness and skin depth 

MAX Number of sampling points 

ITER Running count of the number of iterations 

CYCLES Maximum number of iterations allowed 
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SINHH , COSHH 

FNCTZ 

XFORM 

TRCOPR_ 

FFT 


Subroutines and Functions 

Complex hyperbolic sin and cosine for 
a complex argument 

The residual vector 

Original transformation for the electric 
field 

Subroutines dependent upon cell geometry 
Fast fourier transform 
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START fnct:-: 



STOP 


Figure 17. Flow Chart 

























C **************** SECANT32.FOR **************************** 

C 

C SECANT METHOD APPLIED TO ASSURE CONVERGENCE 
C 

C OCTOBER 23,1983 BY ROBERT MIDDELVEEN 
C 

C DIMENSION ALL ARRAYS 

COMPLEX E(32) , FI { 32 ) , FIM1 ( 32 ) ,JC(32),G(32) 

COMPLEX RI,CREF,HI,EHOLD(32) 

COMPLEX GUESS( 32 ) 

COMPLEX J,CK1(32) ,CK2(32) ,CK(32) , Z , T1 , SINHH , COSHH 

REAL K,K2, RSI, RS2,SKIND1, RATIO 

INTEGER ITER, CYCLES 

LOGICAL CONVERGED 

COMMON RI,CREF,HI ,EINC 

COMMON JC , G , J , Z , CK 

COMMON N,Nl,IW,MAX,W, UU , STH , DR, REF, B, ITER 
C A = FLOQUET CELL DIMENSION 
C B = STRIP SIZE 

WRITE(*,*) ' ** ENTERING MAIN PROGRAM ** ' 

WRITE ( * , * ) ' HOW MANY ITERATIONS DO YOU WISH TO PERFORM?' 

READ (*,*) CYCLES 

WRITE ( * , * ) ' INPUT FLOQUET CELL SIZE, STRIP SIZE ' 

WRITE ( * , * ) ' NORMALIZED IN WAVELENGTHS ' 

READ ( * , * ) A , B 
C FREQ = FREQUENCY IN HZ 

WRITE( * , * ) ' INPUT FREQUENCY IN HZ' 

READ ( * , * ) FREQ 

C MAX = FFT SIZE = NUMBER OF SAMPLES PER CELL 
C IW = LOG2 ( MAX ) ; i.e. MAX = 2**IW 
MAX=32 
IW= 5 

C CIG = CONDUCTANCE OF THE STRIP 

WRITE ( * , * ) ' INPUT CONDUCTANCE OF COATING IN SIEMANS ' 

READ( * , * ) SIG1 

WRITE ( * , * ) ' INPUT CONDUCTANCE OF SUBSTRATE IN SIEMANS' 

READ ( * , * ) SIG2 

WRITE( * , * ) ' INPUT THE THICKNESS OF THE COATING IN METERS' 
READ( * , * ) THICKNESS 
C TH = THETA ANGLE OF INCIDENCE 
C PH = PHI ANGLE OF INCIDENCE 

WRITE ( * , * ) ' INPUT THETA ANGLE, PHI ANGLE' 

487 READ ( * , * ) TH , PH 
C INITIALIZE ROUTINE CONSTANTS 
PI=3. 141593 
TPI=2*P I 
C=2 . 997956E8 
UU=4. 0E-7*PI 
EP=8.854E-12 
ETA=SQRT{ UU/EP) 

ITER=0 
J=(0. 0,1.0) 

ALAMB=C/FREQ 
RD=180. 0/PI 
DR=1 ,0/RD 
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CONVERGED=. FALSE. 

C CALCULATE THE INCIDENT ELECTRIC AND MAGNETIC FIELD COMPONENTS 
C FOR THE ELECTRIC FIELD PARALLEL TO THE WIRES, i.e. NO CROSS- 
C POLARIZATION INCLUDED 

IF (PH. LT. 45.0) EINC =1.0 
IF (PH. GE. 45.0) EINC=COS(TH*DR) 

IF (PH. LT. 45.0) STH=0 . 0 

IF (PH. GT. 45.0) STH=TH 

H 1=1. 0/ETA* ( COS (TH*DR)+J*SIN(TH*DR) ) 

IF(PH.LT.45.0) HI=1.0/ETA*C0S(TH*DR) 

C 

C THIS SECTION COMPUTES THE INTERNAL IMPEDANCE FOR A COATED WIRE 

C 

Tl=(l. 0,1.0) *SQRT (PI*FREQ*UU*SIG1 ) 

RSI =SQRT (PI*FREQ*UU/SIG1 ) 

RS2=SQRT(PI*FREQ*UU/SIG2 ) 

SKIND1 = 1/ SQRT ( PI * FREQ*UU*S IG1 ) 

RATIO=TH ICKNESS/ S KINDI 
IF ( RATIO. GT. 4. 00) THEN 

Z=SQRT ( TPI*FREQ*UU/2 .0/SIG1 ) * ( 1 .0, 1.0)/TPI/B 
ELSE 

Z=(1.0,1.0)*( SINHH ( T1*TH ICKNESS )+(RS2/RSl) *COSHH ( T1*THICKNESS ) ) 
Z=Z/ (C0SHH(T1*THICKNESS)+(RS2/RS1)*SINHH(T1*THICKNESS) ) 
Z=Z*RSl/TPI/B 
END IF 

WRITE (*,15) RATIO, Z 

15 FORMAT ('-',' RATIO = ',F10.5,' Z = ( ',E10.2,' , ',E10.2,' )') 

C 

C 

C 

C CALCULATE THE NUMBER OF SAMPLES ON THE STRIP AND IN THE APERTURE 
TAU=A-B 

N=IFIX ( TAU/A*FLOAT ( MAX ) ) 

WRITE (*,*) ' N= ' , N 

Nl=N+l 

WRITE (*,30) A,B,TAU,FREQ,TH,N,MAX 

30 FORMAT ( ' — ' ,1E10.5, 1 ',1E10.5,' ',1E10.5, ’ ' , E10. 3 , 1E10 . 4 , 21 10 ) 

IF ( N1 . GT . MAX ) GOTO 998 

K=TPI/ALAMB 

K2=K*K 

SK=K*SIN( TH*DR ) *COS ( PH*DR ) 

SSK=K*SIN(TH*DR)*SIN(PH*DR) 

W=TPI*FREQ 

C CALCULATE GREEN FUNCTION TRANSFORM 
DO 40 1=1, MAX 
IF(I.GT. MAX/2+1) GOTO 50 
U=TPI*(I-1)/A-SK 
GOTO 60 

50 U=TPI* ( I-MAX-1 )/A-SK 

60 U=U*U+SSK*SSK 

IF ( U . GE . K2 ) GOTO 70 
G ( I ) =-J*SQRT( K2-U ) 

GOTO 44 

70 G ( I ) =-SQRT ( U-K2 ) 

44 G(I)=G(I)-SSK*SSK/G(I) 



40 CONTINUE 
C INITIAL E FIELD ESIMATE 
DO 320 1=1/ MAX 

320 E(I)=(1. 0,0.0) 

CALL TRCOPR ( E , N 1 ) 

DO 321 1=1, MAX 

321 GUESS(I)=E(I)+(0. 1,0.0) 

C NOTE ITERATIVE FORM USED IN THIS PROGRAM IS X=AX+B 

C CALCULATE B PORTION OF ITERATIVE EQUATION 
DO 110 1=1, MAX 
110 CK1(I)=HI*J*W*UU 
CALL FFT ( CK1 , IW ) 

DO 120 I =N 1 , MAX 
120 CK2( I )=HI*W*UU/J 
DO 140 1=1, N 
140 CK2(I)=(0. 0,0.0) 

CALL FFT ( CK2 , IW ) 

DO 130 1=1, MAX 

130 CK(I)=(CK1(I)+CK2(I) )/G(l) 

C 

C THE FOLLOWING SECTION IMPLEMENTS THE SECANT METHOD 

C 

WRITE(*,*) 'SECANT ALGORITHM APPLIED' 

81 CONTINUE 

ITER=ITER+1 
DO 357 1=1, MAX 
F I ( I ) =E ( I ) 

FIM1 ( I ) =GUESS ( I ) 

EHOLD ( I ) =E ( I ) 

357 CONTINUE 

CALL FNCTZ(FI , MAX, CONVERGED) 

IF (CONVERGED) GOTO 259 

CALL FNCTZ(FIM1 , MAX, CONVERGED) 

IF (ITER. GT. CYCLES. OR. CONVERGED) GOTO 259 
DO 358 1=1, MAX 

E(I)=E(I)-FI(I)*( ( E ( I ) -GUESS ( I ) ) / ( FI ( I ) -FIM1 ( I ) )) 

358 CONTINUE 

DO 359 1=1, MAX 
GUESS ( I ) =EHOLD ( I ) 

359 CONTINUE 
GOTO 81 

C 

C 

C 

C THE SECANT METHOD ENDS HERE 

C 
C 
C 

259 CONTINUE 

C 

C 

WRITE ( * , * ) ' HIT 'RETURN" TO CONTINUE’ 

READ ( * , * ) 

OPEN ( 10 , FI LE= ' SEC320UT ' , STATUS= ' NEW ' ) 

WRITE( * , * ) ' ELECTRIC FIELD MAGNITUDE' 
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WRITE(10,*) ' ELECTRIC FIELD MAGNITUDE' 

DO 261 1=1, MAX 

BINDEX = FLOAT( I -1 ) /FLOAT ( MAX-1 ) *A 
WRITE(10,*) ' ', BINDEX, 1 ',CABS(E(I)) 

WRITE( * , * ) ' ' , ' BINDEX = ’, BINDEX, 1 ! E 1 = ',CABS(E(I)) 

261 CONTINUE 

WRITE ( * , * ) 1 HIT 'RETURN" TO CONTINUE 1 
READ (*,*) 

WRITE ( * , * ) ’ STRIP CURRENT 1 

WRITE(10,*) 1 STRIP CURRENT* 

DO 262 I=N1 , MAX 

BINDEX=FLOAT( 1-1 ) /FLOAT (MAX-1 )*A 
WRITEdO,*) 1 BINDEX, 1 1 , CABS ( JC ( I ) ) 

WRITE(*,*)’ ’, ’BINDEX = ’, BINDEX, 1 !JC! = 1 , CABS ( JC ( I ) ) 

262 CONTINUE 
C 

C 

WRITE ( * , 2 60 ) 

260 FORMAT ( 1 - 1 , 10X, 1 TIMELY EXIT 1 ) 

GOTO 9999 
998 WRITE( * , 99 ) 

99 FORMAT ( 1 - 1 , 1 ERROR IN N 1 ) 

9999 STOP 
END 
C 
C 
C 

C THIS SUBROUTINE PRODUCES THE VECTOR WE WANT TO ZERO 

C 
C 
C 

SUBROUTINE FNCTZ ( E , MAX , CONVERGED ) 

COMPLEX E( 32 ) , HOLD( 32 ) 

LOGICAL CONVERGED 
DO 1234 1=1, MAX 
HOLD ( I ) =— E ( I ) 

1234 CONTINUE 

CALL XFORM(E, CONVERGED) 

DO 1235 1=1, MAX 
E ( I ) =E ( I ) +HOLD ( I ) 

1235 CONTINUE 
RETURN 
END 

C 

C 

C 

C THIS SUBROUTINE IMPLEMENTS THE TRANSFORMATION ON THE E FIELD 

C 
C 
C 

SUBROUTINE XFORM ( E , CONVERGED ) 

COMPLEX E(32) , G ( 32 ) , JC ( 3 2 ) , CK ( 32 ) 

COMPLEX RI , CREF , HI , J , Z 
REAL REF,REFM1,REFM2 
LOGICAL CONVERGED 
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COMMON RI/CREF/HI/EINC 
COMMON JC,G,J,Z,CK 

COMMON N,N1,IW,MAX,W,UU,STH,DR,REF,B,ITER 
C CALCULATE FIELD ON STRIP DUE TO FINITE CONDUCTIVITY 

CALL TRCOPR3(E,JC,Nl,MAX,Z,B) 

C START BY PERFORMING THE INITIAL TRANSFORMATION 

CALL FFT ( E / IW ) 

DO 100 1=1, MAX 
100 E ( I ) =CONJG( E ( I ) *G ( I ) ) 

C PERFORM INVERSE TRANSFORM OF ( G* E ) 

CALL FFT ( E , IW ) 

C PERFORM THE TRUNCATION OPERATION T ( G* E ) 

CALL TRCOPRC ( E , N ) 

CALL TRC0PR4 (E,JC,N1 ,MAX,J ,W,UU,HI ) 

C PERFORM INVERSE TRANSFORMATION ON T ( G* E ) 

CALL FFT ( E , I W ) 

C PERFORM T(G*E)/G AND ADD CONSTANT "B" 

DO 170 1=1, MAX 
E(I)=E(I)/G(I) +CK ( I ) 

170 CONTINUE 
C 

C 

C 

C CALCULATE REFLECTION COEFFICIENT 
TOL=0 .0001 
REFM2=REFMl 
REFM1=REF 

RI=J*SIN(STH*DR) /COS ( STH*DR ) 

CREF=(E(1)/MAX+EINC)+J*SIN(STH*DR)*ABS(1.0— ABS(E(1)/MAX+EINC) ) 
CREF=CREF/ ( COS ( STH*DR ) +J*S IN ( STH*DR ) ) 

REF=CABS ( CREF ) 

IF (ABS(REF-REFM1 ) . LT . TOL . AND . ABS ( REF-REFM2 ) .LT.TOL) THEN 
WRITE ( * , * ) ' ITER= ’ , ITER , ' REF= ' , REF , ' CREF= ' , CREF 
CONVERGED= . TRUE . 

END IF 
C 
C 

DO 171 1=1, MAX 
E ( I ) =CON JG ( E { I ) ) /MAX 

171 CONTINUE 

C PERFORM INVERSE TRANSFORMATION TO OBTAIN FIRST ITERATED 
C ELECTRIC FIELD 

CALL FFT ( E , IW ) 

DO 200 1=1, MAX 
200 E ( I ) =CONJG ( E ( I ) ) 

RETURN 

END 

C 

C 

C 

SUBROUTINE TRCOPR(E,N) 

COMPLEX E( 32 ) 

DO 400 I=N , 32 
E(I )=(0. 0,0.0) 

400 CONTINUE 



RETURN 

END 

C 

C 

c 

SUBROUTINE TRCOPRC(E,N) 

COMPLEX E( 32 ) 

DO 401 1=1, N 
E(I)=(0. 0,0.0) 

401 CONTINUE 
RETURN 
END 

C 

C 

C 

SUBROUTINE TRCOPR3 ( E , JC , N1 , MAX , Z , B) 

COMPLEX E ( 32 ) ,JC(32) ,Z 
DO 402 I=N1 , MAX 
E ( I ) =-J C ( I ) *Z *B 

402 CONTINUE 
RETURN 
END 

C 

C 

C 

SUBROUTINE TRCOPR4 ( E , JC , N1 , MAX , J , W , UU , HI ) 

COMPLEX E(32),JC(32),J,HI 
DO 403 I =N 1 , MAX 
E ( I ) =CONJG ( E ( I ) ) /MAX 

C CALCULATE THE CURRENT DENSITY ON THE STRIP 
JC( I ) -E ( I ) * J/W/UU— HI 

403 CONTINUE 
RETURN 
END 

C 

C 

C 

SUBROUTINE FFT ( A , M ) 

C THIS IS THE FFT SUBROUTINE CALLED FOR FROM THE MAIN PROGRAM 
COMPLEX A( 32) ,U,W,T 
N=2**M 
NV2=N/2 
NMl =N— 1 
J = 1 

DO 7 1=1, NMl 
IF(I.GE.J) GOTO 5 
T=A( J) 

A ( J ) =A ( I ) 

A ( I ) =T 

5 K=NV2 

6 IF(K.GE.J) GOTO 7 
J=J-K 
K=K/2 
GOTO 6 
J=J+K 


7 



o n n n o o 


PI=3. 141592653589793 

DO 20 L=1 / M 

LE=2**L 

LEl=LE/2 

U=( 1.0, 0.0) 

W=CMPLX(C0S(PI/LE1 ) ,SIN{PI/LE1)) 
DO 20 J=1,LE1 
DO 10 I=J,N,LE 
IP= I+LE1 
T=A ( IP ) *13 
A ( I P ) =A ( I ) -T 
10 A ( I ) =A ( I ) +T 

20 U=U*W 

RETURN 
END 


COMPLEX FUNCTION SINHH(X) 
COMPLEX X 

SINHH=0.5*(CEXP(X)-CEXP{ -X ) ) 
END 


COMPLEX FUNCTION COSHH(X) 
COMPLEX X 

COSHH=0 . 5* ( CEXP ( X ) +CEXP ( -X ) ) 
END 



APPENDIX B 


The purpose of this appendix is to solve for the 
magnetic field intensity as a function of magnetic vector 
potential/ this derivation is included in reference 7. Given 
equation 1/ it is desired to derive equation 4. For this 
derivation all sources will be considered sinusoidal, E and 
IT will be phasor quantities and the sinusoidal steady state 
versions of Maxwell's equations will be used along with 
three vector identities. Considering X and B as arbitrary 
vectors and C as an arbitrary scalar: 


7 X 

>| 

+ 

ro| 

) = 7 x X + 7 x b 

[56] 

7 x 

7 C = 

0 

[57] 

7 x 

V x “ 

_ 2_ 

= V ( F-A) - 7a 

[58] 

Starting with 

the relation of electric and magnetic 

fields 


given by maxwell's equation: 

7 X TT = ju>€ ¥ [59] 

Equation 1, repeated below, 

e" = -1/C ( 7 X T ) [1] 

is substituted in equation 59 becoming: 

PXH - + ju( F X 7 ) = 0 [60] 

Making use of identity 56 this can be written as: 

7 X ( "H + j w “ ) = 0 [61] 
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Making use of identity 57 this can be written again as: 


V X ( H + j(j F ) = 

v x ( -vK ) 



[62] 

The magnetic scalar 

potential is denoted 

by 

0m 

in this 

equation. From 62 the 

following implication 

can 

be 

made : 

TT + j t*>F = -V<Pm 




[63] 


To specify H completely from F it is necessary to find 
the relation between 0m and F. Taking the curl of equation 1 
results in: 

V X E = -1/e V X P X P [64] 

Making use of identity 58 this can be written as: 

V x E = -1/e [? ( 7-F) - V 2 ~ ] [65] 

To specify any vector both its curl and divergence must be 
known. The curl of F is already defined in equation 1. To 
specify the divergence of F, the Lorentz gauge condition is 
applied so that equation 65 is simplified. The divergence of 
F is defined by: 

P*F = -j U€M<t> m [66] 

The scalar magnetic potential is now given in terms of 
vector magnetic potential by: 

<t> m = -1/j <*>»e ( V-F ) [67] 

This value can be substituted back into equat'ion 62 
resulting in eqation 4: 

F = -j wF + l/jc one 7( P*F ) 


[ 4 ] 


APPENDIX C 


The following proof is also found in the reference of 
Brand [6], It is a proof that the Newton-Raphson method is 
identical to the contraction corrector method. The starting 
point is the definition of the Newton-Raphson method. 

x i+1 = x i - f(x i )/f'(x i ) [68] 
where it is desired that 


f(x i ) = 0 [69] 

and that equivalently 

g(x i )-x i =0 [70] 

equation 68 can now be manipulated as follows: 




[73] 

[74] 
g(x i ) [75] 
gCx^ [76] 


The optimum correction factor is defined by: 


R = g ' ( x ^ ) / [ g ' (x.^) - 1 ] 


[77] 
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So that equation 75 becomes the definition of the optimum 
contraction corrector when the substitution stated in 
equation 77 is made as follows: 

x i + l= R x i + ( 1 - R ). g(x i ) [78] 



APPENDIX D 


The solution for the internal impedance of the coated 
conductor (refer to reference 2) is found by solving for the 
distribution equations in both media and then matching at 
the boundary between the two (see Figure 19). The solution 
for either material, assuming an electric field with only a 
z component, is of the following form: 

i z = i 0 EXP(-x/6) EXP(-jx/«) [79] 
Figure 19 Coated Conductor 



There can be no positive exponential for the substrate, 

however, because the current becomes zero for large values 

of x. The current density for the substrate then becomes: 

i „ = C EXP ( -T _ x) [80] 

z2 ^ 

where T 2 = ( 1 + j ) = (1 + j ) 7T f t^<f 2 

b 2 
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The current density for the coating has both the positive 
and negative exponetial: 

i zl = D EXPf-Tj x) + E EXPC+Tj^ x) [81] 

where Tj = ( 1 + j ) = (1 + j 

6 

It is more convenient to express this expression in terms of 

the equivalent hyperbolic functions. Then i , becomes: 

zl 

i zl = A SINH(T 1 x) + B COSHCTj x) [82] 

The electric and magnetic fields now need to be matched 
at the boundary. These fields can be found using the 


following relations: 

E z = i z /a [83] 

Hy = 1 d(E^) = <rd(E^) [84] 

j U)H dx T 2 dx 

Solving for the electric and magnetic fields yields: 

e z2 = £ EXP(-T 2 x) [85] 

a 2 

E zl = £( A SINH(T 1 x) + B COSH ( T^ x)) [86] 

... ^ 

H y2 = -C EXP(-T 2 x) [87] 

t 2 

Hyl = ].( A COSH(T 1 x) + B SINH^ X)) [88] 

T 1 


The tangential electric and magnetic fields are continuous 
across the boundaries/ yielding two equations and three 



unknowns. Combining equations will obtain the ratio/ B/A: 
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E zl = E z2 ' H y 1 = H y2 @ x = d 

B = - [SINH(T 1 d) + (T 2 a /T 1 a )COSH(T-. d)] 

A [ COSH { 2 d) + ( T 2 a £)SINH(tJ d)] 

The total current is obtained from the relationship: 

7 = 7x7 


[89] 

[90] 


[91] 


Solving fo J : 

J z = -H y [92] 

The impedance per square can be obtained from observing the 

ratio of the electric field and the current density at x=0: 

Z = = - E = - B T 1 at x = 0 [93] 

J z H y A ^ 

Using equation 90 and substituting for T1 and T2 yields the 
ratio of the internal impedance to the surface resistivity 
of the coating. 

Z = (1 + j ) [SINH ( T-| d) + (R ,/R ^ )C0SH(T 1 d ) ] [94] 

R sl [COSH(Tj d) + (Rg2/Rsl )SINH(T l 
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